data(beets)
head(beets)
beet0<-lmer(sugpct~block+sow+harvest+(1|block:harvest), data=beets, REML=FALSE)
beet_no.harv <- update(beet0, .~.-harvest)
PBmodcomp(beet0, beet_no.harv, nsim=20)
## Vanilla
PBmodcomp(beet0, beet_no.harv)
## Simulate reference distribution separately:
rr <- PBrefdist(beet0, beet_no.harv, nsim=20)
PBmodcomp(beet0, beet_no.harv, ref=rr)
Run the code above in your browser using DataLab