> library(foreign)
> ds <- read.dta("c:/amenorrhea.dta")
> attach(ds)
> time2 <- time^2
> trt.time <- trt*time
> trt.time2 <- trt*time2
> library(lme4)
> model1 <- glmer(y ~ time + time2 + trt.time + trt.time2 + (1 | id),
+ family=binomial, nAGQ=50, na.action=na.omit)
> summary(model1)
Generalized linear mixed model fit by the adaptive Gaussian Hermite approximation
Formula: y ~ time + time2 + trt.time + trt.time2 + (1 | id)
AIC BIC logLik deviance
4392 4430 -2190 4380
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 4.2461 2.0606
Number of obs: 3616, groups: id, 1151
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.73360 0.26811 -13.926 < 2e-16 ***
time 1.11339 0.25062 4.443 8.89e-06 ***
time2 -0.04234 0.05140 -0.824 0.41005
trt.time 0.54365 0.17721 3.068 0.00216 **
trt.time2 -0.10529 0.04621 -2.278 0.02271 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) time time2 trt.tm
time -0.859
time2 0.759 -0.971
trt.time -0.014 -0.360 0.441
trt.time2 0.015 0.342 -0.459 -0.955
> model2 <- glmer(y ~ time + time2 + (1 | id), family=binomial,
+ nAGQ=50, na.action=na.omit)
> summary(model2)
Generalized linear mixed model fit by the adaptive Gaussian Hermite approximation
Formula: y ~ time + time2 + (1 | id)
AIC BIC logLik deviance
4402 4426 -2197 4394
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 4.2474 2.0609
Number of obs: 3616, groups: id, 1151
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.72846 0.26757 -13.935 < 2e-16 ***
time 1.38741 0.23328 5.947 2.72e-09 ***
time2 -0.09595 0.04555 -2.107 0.0352 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) time
time -0.926
time2 0.863 -0.981
> anova(model1,model2)
Data:
Models:
model2: y ~ time + time2 + (1 | id)
model1: y ~ time + time2 + trt.time + trt.time2 + (1 | id)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
model2 4 4401.7 4426.5 -2196.9
model1 6 4392.4 4429.6 -2190.2 13.319 2 0.001282 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Marginal Logistic Regression Model
> library(geepack)
> summary(geeglm(y ~ time + time2 + trt.time + trt.time2,
+ id = id, waves=time, family=binomial("logit"),
+ corstr="unstructured", std.err="san.se"))
Call:
geeglm(formula = y ~ time + time2 + trt.time + trt.time2, family = binomial("logit"),
id = id, waves = time, corstr = "unstructured", std.err = "san.se")
Coefficients:
Estimate Std.err Wald Pr(>|W|)
(Intercept) -2.24079 0.17657 161.051 < 2e-16 ***
time 0.69794 0.15813 19.482 1.02e-05 ***
time2 -0.03137 0.03176 0.975 0.32342
trt.time 0.33801 0.10967 9.499 0.00206 **
trt.time2 -0.06903 0.02841 5.902 0.01513 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Estimated Scale Parameters:
Estimate Std.err
(Intercept) 1 0.02834
Correlation: Structure = unstructured Link = identity
Estimated Correlation Parameters:
Estimate Std.err
alpha.1:2 0.3438 0.03662
alpha.1:3 0.2594 0.03591
alpha.1:4 0.2730 0.03408
alpha.2:3 0.4260 0.03579
alpha.2:4 0.3855 0.03450
alpha.3:4 0.4993 0.03607
Number of clusters: 1151 Maximum cluster size: 4
Conditional Logistic Regression Model
> clogit(y ~ time + time2 + trt.time + trt.time2 + strata(id))
Call:
clogit(y ~ time + time2 + trt.time + trt.time2 + strata(id))
coef exp(coef) se(coef) z p
time 0.7587 2.14 0.3766 2.014 0.044
time2 0.0367 1.04 0.0742 0.494 0.620
trt.time 1.2107 3.36 0.5324 2.274 0.023
trt.time2 -0.2235 0.80 0.1040 -2.148 0.032
Likelihood ratio test=426 on 4 df, p=0 n= 3616, number of events= 1231
(988 observations deleted due to missingness)
> library(foreign)
> 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")
> attach(dslong)
> time <- as.numeric(visit!=0)
> ltime <- time
> ltime[visit==0] <- log(8)
> ltime[visit!=0] <- log(2)
> library(lme4)
> summary(glmer(y ~ offset(ltime) + time + trt + time:trt + (1 + time | id),
+ family=poisson, nAGQ=50, na.action=na.omit))
Generalized linear mixed model fit by the adaptive Gaussian Hermite approximation
Formula: y ~ offset(ltime) + time + trt + time:trt + (1 + time | id)
AIC BIC logLik deviance
788 813.7 -387 774
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 0.49620 0.70441
time 0.22791 0.47740 0.075
Number of obs: 295, groups: id, 59
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.068437 0.139439 7.662 1.82e-14 ***
time 0.003752 0.107261 0.035 0.9721
trt 0.052449 0.192112 0.273 0.7848
time:trt -0.305700 0.149259 -2.048 0.0405 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) time trt
time -0.049
trt -0.726 0.035
time:trt 0.035 -0.719 -0.044
> summary(glmer(y ~ offset(ltime) + time + trt + time:trt + (1 + time | id),
+ family=poisson, subset=(id!=49), nAGQ=50, na.action=na.omit))
Generalized linear mixed model fit by the adaptive Gaussian Hermite approximation
Formula: y ~ offset(ltime) + time + trt + time:trt + (1 + time | id)
Subset: id != 49
AIC BIC logLik deviance
755.7 781.4 -370.9 741.7
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 0.44888 0.66999
time 0.21497 0.46365 -0.029
Number of obs: 290, groups: id, 58
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.067267 0.133261 8.009 1.16e-15 ***
time 0.011989 0.105026 0.114 0.9091
trt -0.006744 0.185153 -0.036 0.9709
time:trt -0.346202 0.147911 -2.341 0.0193 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) time trt
time -0.140
trt -0.720 0.101
time:trt 0.099 -0.710 -0.136
> library(foreign)
> ds <- read.dta("c:/arthritis.dta")
> dslong <- reshape(ds, idvar="id", varying=c("y1","y2","y3","y4"),
+ v.names="y", timevar="time", time=1:4, direction="long")
> dslong <- dslong[order(dslong$id, dslong$time),]
> dslong <- dslong[complete.cases(dslong$y),]
> attach(dslong)
> month <- 2*(time-1)
> sqrtmonth <- month^0.5
> library(ordinal)
> summary(clmm(ordered(y) ~ trt + sqrtmonth + trt:sqrtmonth,
+ random = factor(id), nAGQ = 50, Hess = TRUE))
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 50 quadrature points
Call:
clmm(location = ordered(y) ~ trt + sqrtmonth + trt:sqrtmonth,
random = factor(id), Hess = TRUE, nAGQ = 50)
Random effects:
Var Std.Dev
factor(id) 4.107657 2.026736
Location coefficients:
Estimate Std. Error z value Pr(>|z|)
trt -0.1270 0.3160 -0.4021 0.6876464
sqrtmonth -0.3759 0.0887 -4.2395 2.2402e-05
trt:sqrtmonth -0.3706 0.1251 -2.9622 0.0030547
No scale coefficients
Threshold coefficients:
Estimate Std. Error z value
1|2 -4.9720 0.2942 -16.8993
2|3 -1.9150 0.2404 -7.9651
3|4 1.0080 0.2345 4.2976
4|5 3.6997 0.2787 13.2730
log-likelihood: -1440.072
AIC: 2896.145
Condition number of Hessian: 141.6181