# NOT RUN {
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, cl=1)
anova(sug, sug.h)
PBmodcomp(sug, sug.s, nsim=50, cl=1)
## 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, cl=1)
anova(sug, sug.s)
PBmodcomp(sug, sug.s, nsim=50, cl=1)
## 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, cl=1)
anova(glm.D93, glm.D93.t, test="Chisq")
PBmodcomp(glm.D93, glm.D93.t, nsim=50, cl=1)
## 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, cl=2)
# }
# 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, cl=1)
## The same test using a restriction matrix
L <- cbind(0,1)
PBmodcomp(fmLarge, L, cl=1)
## Vanilla
PBmodcomp(beet0, beet_no.harv, nsim=1000, cl=1)
## Simulate reference distribution separately:
refdist <- PBrefdist(beet0, beet_no.harv, nsim=1000)
PBmodcomp(beet0, beet_no.harv, ref=refdist, cl=1)
## 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)
# }
# NOT RUN {
# }
Run the code above in your browser using DataCamp Workspace