> library(foreign)
> library(MASS)
> library(lme4)
> ds <- read.dta("c:/epilepsy.dta")
> dslong <- reshape(ds, idvar="id", varying=c("y0","y1","y2","y3","y4"),
+ v.names="y", timevar="visit", time=0:4, direction="long")
> dslong <- subset(dslong,id!=49)
> attach(dslong)
> time <- as.numeric(visit!=0)
> ltime <- time
> ltime[visit==0] <- log(8)
> ltime[visit!=0] <- log(2)
> ## PQL estimation
> summary(glmmPQL(y ~ offset(ltime) + time + trt + time:trt,
+ random = ~1 + time | id, family=poisson))
Loading required package: nlme
Attaching package: 'nlme'
The following object(s) are masked from 'package:lme4':
BIC, fixef, lmList, ranef, VarCorr
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
Linear mixed-effects model fit by maximum likelihood
Data: NULL
AIC BIC logLik
NA NA NA
Random effects:
Formula: ~1 + time | id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.6348082 (Intr)
time 0.3826100 0.164
Residual 1.3920269
Variance function:
Structure: fixed weights
Formula: ~invwt
Fixed effects: y ~ offset(ltime) + time + trt + time:trt
Value Std.Error DF t-value p-value
(Intercept) 1.1212102 0.1326893 230 8.449889 0.0000
time -0.0028608 0.1049100 230 -0.027269 0.9783
trt -0.0225051 0.1845751 56 -0.121929 0.9034
time:trt -0.3005242 0.1491862 230 -2.014423 0.0451
Correlation:
(Intr) time trt
time -0.092
trt -0.719 0.066
time:trt 0.065 -0.703 -0.091
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.1218196 -0.5883254 -0.1072718 0.3428562 5.2983667
Number of Observations: 290
Number of Groups: 58
> ## Laplace approximation
> summary(glmer(y ~ offset(ltime) + time + trt + time:trt + (1 + time | id),
+ family=poisson, nAGQ=1))
Generalized linear mixed model fit by the Laplace approximation
Formula: y ~ offset(ltime) + time + trt + time:trt + (1 + time | id)
AIC BIC logLik deviance
771.4 797 -378.7 757.4
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 0.45182 0.67218
time 0.21489 0.46356 0.049
Number of obs: 290, groups: id, 58
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.069379 0.133603 8.004 1.20e-15 ***
time 0.007766 0.104959 0.074 0.9410
trt -0.008019 0.185643 -0.043 0.9655
time:trt -0.345976 0.147852 -2.340 0.0193 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) time trt
time -0.076
trt -0.720 0.055
time:trt 0.054 -0.710 -0.073
> library(foreign)
> library(MASS)
> library(lme4)
> ds <- read.dta("c:/amenorrhea.dta")
> attach(ds)
> time2 <- time^2
> trt.time <- trt*time
> trt.time2 <- trt*time2
> ## PQL estimation
> summary(glmmPQL(y ~ time + time2 + trt.time + trt.time2,
+ random = ~ 1 | id, family = binomial, na.action=na.omit))
Loading required package: nlme
Attaching package: 'nlme'
The following object(s) are masked from 'package:lme4':
BIC, fixef, lmList, ranef, VarCorr
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
Linear mixed-effects model fit by maximum likelihood
Data: NULL
AIC BIC logLik
NA NA NA
Random effects:
Formula: ~1 | id
(Intercept) Residual
StdDev: 2.048812 0.7574991
Variance function:
Structure: fixed weights
Formula: ~invwt
Fixed effects: y ~ time + time2 + trt.time + trt.time2
Value Std.Error DF t-value p-value
(Intercept) -3.385835 0.21018273 2461 -16.109008 0.0000
time 1.035071 0.19872753 2461 5.208495 0.0000
time2 -0.034102 0.04085522 2461 -0.834713 0.4040
trt.time 0.508912 0.15324846 2461 3.320829 0.0009
trt.time2 -0.098122 0.03831391 2461 -2.560994 0.0105
Correlation:
(Intr) time time2 trt.tm
time -0.831
time2 0.737 -0.971
trt.time -0.008 -0.394 0.459
trt.time2 0.010 0.375 -0.480 -0.952
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.9779646 -0.5580611 -0.2761477 0.5411420 4.9843269
Number of Observations: 3616
Number of Groups: 1151
> ## Laplace approximation
> summary(glmer(y ~ time + time2 + trt.time + trt.time2 + (1 | id),
+ family = binomial, nAGQ=1, na.action=na.omit))
Generalized linear mixed model fit by the Laplace approximation
Formula: y ~ time + time2 + trt.time + trt.time2 + (1 | id)
AIC BIC logLik deviance
3927 3964 -1957 3915
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 4.3366 2.0825
Number of obs: 3616, groups: id, 1151
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.76280 0.26910 -13.983 < 2e-16 ***
time 1.12165 0.25150 4.460 8.2e-06 ***
time2 -0.04331 0.05158 -0.840 0.40109
trt.time 0.54717 0.17829 3.069 0.00215 **
trt.time2 -0.10565 0.04643 -2.275 0.02289 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) time time2 trt.tm
time -0.858
time2 0.759 -0.971
trt.time -0.014 -0.361 0.442
trt.time2 0.015 0.343 -0.460 -0.955