### R code from vignette source 'glmm2.Rnw'

###################################################
### code chunk number 1: zero
###################################################
options(show.signif.stars=FALSE,width=80)
set.seed(10131985)


###################################################
### code chunk number 2: two
###################################################
url <-"http://tinyurl.com/3dmnytg"
data <- read.csv(url)
data$Nchild <- factor(data$Nchild)
data$Area <-factor(data$Area)
data$Cage <- data$Age - mean(data$Age)
require(car)
library(effects)
require(lme4)
print(
 xyplot(Cuse ~ Age|Nchild*Urban, data, type=c("g", "r", "smooth", "p"))
)


###################################################
### code chunk number 3: three
###################################################
m1 <- glm(Cuse ~ (Nchild + poly(Cage, 2, raw=TRUE) + Urban)^2, data=data, family=binomial)
Anova(m1)
Anova(m2 <- update(m1, ~ Nchild*poly(Cage, 2, raw=TRUE) + Urban))
print(plot(effect("Nchild*poly(Cage, 2, raw=TRUE)", m2), multiline=TRUE, type="link"))


###################################################
### code chunk number 4: threea
###################################################
data$kids <- factor(ifelse(data$Nchild == 0, "n", "y"))
m3 <- update(m2, ~ kids*poly(Cage, 2, raw=TRUE) + Urban)
print(plot(effect("kids*poly(Cage, 2, raw=TRUE)", m3), multiline=TRUE, type="link"))
Anova(m3)


###################################################
### code chunk number 5: four
###################################################
g2 <- lmer(Cuse ~ kids*poly(Cage, 2, raw=TRUE) + Urban + (1|Area), data=data, family=binomial)
Anova(g2)
print(plot(effect("kids:poly(Cage, 2, raw=TRUE)", g2), multiline=TRUE, type="link"))


###################################################
### code chunk number 6: glmm2.Rnw:81-83
###################################################
compareCoefs(m3, g2)
exp(fixef(g2))


###################################################
### code chunk number 7: five
###################################################
str(ranef(g2))
par(mfrow=c(1, 2))
plot(density(ff <- ranef(g2)$Area[, 1]), main="Density of Area effects",
   ylim=c(0, 1.2), lwd=2)
xx <- seq(-3.5, 3.5, length=200)
lines(mean(ff) + sd(ff)*xx, dnorm(xx)/sd(ff))
qqnorm(ff, main="QQ-plot of random effects")
qqline(ff)


###################################################
### code chunk number 8: glmm2.Rnw:97-99
###################################################
print(g3 <- update(g2, ~ . - (1|Area) + (Urban -1 | Area)), corr=FALSE)
anova(g2, g3, text="Chisq")