data(beets, package="pbkrtest")
head(beets)
## Linear mixed effects model:
sug <- lmer(sugpct ~ block + sow + harvest + (1|block:harvest), data=beets, REML=FALSE)
sug.h <- update(sug, .~. -harvest)
sug.s <- update(sug, .~. -sow)
anova(sug, sug.h)
PBmodcomp(sug, sug.h, nsim=50)
anova(sug, sug.h)
PBmodcomp(sug, sug.s, nsim=50)
## Linear normal model:
sug <- lm(sugpct ~ block + sow + harvest, data=beets)
sug.h <- update(sug, .~. -harvest)
sug.s <- update(sug, .~. -sow)
anova(sug, sug.h)
PBmodcomp(sug, sug.h, nsim=50)
anova(sug, sug.s)
PBmodcomp(sug, sug.s, nsim=50)
## Generalized linear model
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
d.AD <- data.frame(treatment, outcome, counts)
head(d.AD)
glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())
glm.D93.o <- update(glm.D93, .~. -outcome)
glm.D93.t <- update(glm.D93, .~. -treatment)
anova(glm.D93, glm.D93.o, test="Chisq")
PBmodcomp(glm.D93, glm.D93.o, nsim=50)
anova(glm.D93, glm.D93.t, test="Chisq")
PBmodcomp(glm.D93, glm.D93.t, nsim=50)
## Generalized linear mixed model (it takes a while to fit these)
## Not run:
# (gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
# data = cbpp, family = binomial))
# (gm2 <- update(gm1, .~.-period))
# anova(gm1, gm2)
# PBmodcomp(gm1, gm2)
# ## End(Not run)
## Not run:
# (fmLarge <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
# ## removing Days
# (fmSmall <- lmer(Reaction ~ 1 + (Days|Subject), sleepstudy))
# anova(fmLarge, fmSmall)
# PBmodcomp(fmLarge, fmSmall)
#
# ## The same test using a restriction matrix
# L<-cbind(0,1)
# PBmodcomp(fmLarge, L)
#
# ## Vanilla
# PBmodcomp(beet0, beet_no.harv, nsim=1000)
#
# ## Simulate reference distribution separately:
# refdist <- PBrefdist(beet0, beet_no.harv, nsim=1000)
# PBmodcomp(beet0, beet_no.harv, ref=refdist)
#
# ## Do computations with multiple processors:
# ## Number of cores:
# (nc <- detectCores())
# ## Create clusters
# cl <- makeCluster(rep("localhost", nc))
#
# ## Then do:
# PBmodcomp(beet0, beet_no.harv, cl=cl)
#
# ## Or in two steps:
# refdist <- PBrefdist(beet0, beet_no.harv, nsim=1000, cl=cl)
# PBmodcomp(beet0, beet_no.harv, ref=refdist)
#
# ## It is recommended to stop the clusters before quitting R:
# stopCluster(cl)
# ## End(Not run)
Run the code above in your browser using DataLab