epil.RdThall and Vail (1990) give a data set on two-week seizure counts for 59 epileptics. The number of seizures was recorded for a baseline period of 8 weeks, and then patients were randomly assigned to a treatment group or a control group. Counts were then recorded for four successive two-week periods. The subject's age is the only covariate.
epilThis data frame has 236 rows and the following 9 columns:
ythe count for the 2-week period.
trttreatment, "placebo" or "progabide".
basethe counts in the baseline 8-week period.
agesubject's age, in years.
V40/1 indicator variable of period 4.
subjectsubject number, 1 to 59.
periodperiod, 1 to 4.
lbaselog-counts for the baseline period, centred to have zero mean.
lagelog-ages, centred to have zero mean.
The value of y in row 31 was corrected from 21 to
23 in version 7.3-65.
Thall, P. F. and Vail, S. C. (1990) Some covariance models for longitudinal count data with over-dispersion. Biometrics 46, 657–671.
Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth Edition. Springer.
summary(glm(y ~ lbase*trt + lage + V4, family = poisson,
data = epil), correlation = FALSE)
#>
#> Call:
#> glm(formula = y ~ lbase * trt + lage + V4, family = poisson,
#> data = epil)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.89915 0.04258 44.602 < 2e-16 ***
#> lbase 0.94952 0.04356 21.797 < 2e-16 ***
#> trtprogabide -0.34713 0.06098 -5.693 1.25e-08 ***
#> lage 0.89705 0.11644 7.704 1.32e-14 ***
#> V4 -0.16109 0.05458 -2.952 0.00316 **
#> lbase:trtprogabide 0.56223 0.06350 8.855 < 2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> (Dispersion parameter for poisson family taken to be 1)
#>
#> Null deviance: 2521.75 on 235 degrees of freedom
#> Residual deviance: 869.32 on 230 degrees of freedom
#> AIC: 1647.3
#>
#> Number of Fisher Scoring iterations: 5
#>
epil2 <- epil[epil$period == 1, ]
epil2["period"] <- rep(0, 59); epil2["y"] <- epil2["base"]
epil["time"] <- 1; epil2["time"] <- 4
epil2 <- rbind(epil, epil2)
epil2$pred <- unclass(epil2$trt) * (epil2$period > 0)
epil2$subject <- factor(epil2$subject)
epil3 <- aggregate(epil2, list(epil2$subject, epil2$period > 0),
function(x) if(is.numeric(x)) sum(x) else x[1])
epil3$pred <- factor(epil3$pred,
labels = c("base", "placebo", "drug"))
contrasts(epil3$pred) <- structure(contr.sdif(3),
dimnames = list(NULL, c("placebo-base", "drug-placebo")))
## IGNORE_RDIFF_BEGIN
summary(glm(y ~ pred + factor(subject) + offset(log(time)),
family = poisson, data = epil3), correlation = FALSE)
#>
#> Call:
#> glm(formula = y ~ pred + factor(subject) + offset(log(time)),
#> family = poisson, data = epil3)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.122e+00 2.008e-01 5.587 2.30e-08 ***
#> predplacebo-base 1.108e-01 4.689e-02 2.363 0.018127 *
#> preddrug-placebo -1.037e-01 6.505e-02 -1.594 0.110992
#> factor(subject)2 7.448e-16 2.828e-01 0.000 1.000000
#> factor(subject)3 -3.857e-01 3.144e-01 -1.227 0.219894
#> factor(subject)4 -1.744e-01 2.960e-01 -0.589 0.555847
#> factor(subject)5 1.577e+00 2.197e-01 7.178 7.08e-13 ***
#> factor(subject)6 6.729e-01 2.458e-01 2.738 0.006182 **
#> factor(subject)7 -4.082e-02 2.858e-01 -0.143 0.886411
#> factor(subject)8 1.772e+00 2.163e-01 8.189 2.64e-16 ***
#> factor(subject)9 5.878e-01 2.494e-01 2.356 0.018454 *
#> factor(subject)10 5.423e-01 2.515e-01 2.156 0.031060 *
#> factor(subject)11 1.552e+00 2.202e-01 7.048 1.81e-12 ***
#> factor(subject)12 9.243e-01 2.364e-01 3.910 9.22e-05 ***
#> factor(subject)13 3.075e-01 2.635e-01 1.167 0.243171
#> factor(subject)14 1.212e+00 2.278e-01 5.320 1.04e-07 ***
#> factor(subject)15 1.765e+00 2.164e-01 8.153 3.54e-16 ***
#> factor(subject)16 9.708e-01 2.348e-01 4.134 3.57e-05 ***
#> factor(subject)17 -4.082e-02 2.858e-01 -0.143 0.886411
#> factor(subject)18 2.236e+00 2.104e-01 10.629 < 2e-16 ***
#> factor(subject)19 2.776e-01 2.651e-01 1.047 0.295060
#> factor(subject)20 3.646e-01 2.603e-01 1.401 0.161324
#> factor(subject)21 3.922e-02 2.801e-01 0.140 0.888645
#> factor(subject)22 -8.338e-02 2.889e-01 -0.289 0.772894
#> factor(subject)23 1.823e-01 2.708e-01 0.673 0.500777
#> factor(subject)24 8.416e-01 2.393e-01 3.517 0.000436 ***
#> factor(subject)25 2.069e+00 2.123e-01 9.750 < 2e-16 ***
#> factor(subject)26 -5.108e-01 3.266e-01 -1.564 0.117799
#> factor(subject)27 -2.231e-01 3.000e-01 -0.744 0.456990
#> factor(subject)28 1.386e+00 2.236e-01 6.200 5.66e-10 ***
#> factor(subject)29 1.605e+00 2.227e-01 7.208 5.70e-13 ***
#> factor(subject)30 1.024e+00 2.372e-01 4.317 1.58e-05 ***
#> factor(subject)31 9.259e-02 2.821e-01 0.328 0.742760
#> factor(subject)32 -3.001e-02 2.909e-01 -0.103 0.917814
#> factor(subject)33 4.721e-01 2.597e-01 1.818 0.069089 .
#> factor(subject)34 3.898e-01 2.640e-01 1.477 0.139764
#> factor(subject)35 1.488e+00 2.251e-01 6.614 3.74e-11 ***
#> factor(subject)36 3.609e-01 2.656e-01 1.359 0.174241
#> factor(subject)37 -1.210e-01 2.979e-01 -0.406 0.684646
#> factor(subject)38 1.345e+00 2.283e-01 5.893 3.78e-09 ***
#> factor(subject)39 1.083e+00 2.354e-01 4.601 4.21e-06 ***
#> factor(subject)40 -7.676e-01 3.634e-01 -2.113 0.034642 *
#> factor(subject)41 1.667e-01 2.772e-01 0.601 0.547598
#> factor(subject)42 5.337e-02 2.848e-01 0.187 0.851369
#> factor(subject)43 1.544e+00 2.239e-01 6.896 5.35e-12 ***
#> factor(subject)44 9.616e-01 2.393e-01 4.019 5.85e-05 ***
#> factor(subject)45 1.178e+00 2.326e-01 5.065 4.08e-07 ***
#> factor(subject)46 -5.265e-01 3.355e-01 -1.569 0.116600
#> factor(subject)47 1.054e+00 2.363e-01 4.461 8.17e-06 ***
#> factor(subject)48 -5.265e-01 3.355e-01 -1.569 0.116600
#> factor(subject)49 2.950e+00 2.082e-01 14.173 < 2e-16 ***
#> factor(subject)50 3.898e-01 2.640e-01 1.477 0.139764
#> factor(subject)51 1.039e+00 2.367e-01 4.389 1.14e-05 ***
#> factor(subject)52 5.722e-01 2.548e-01 2.245 0.024746 *
#> factor(subject)53 1.671e+00 2.215e-01 7.543 4.59e-14 ***
#> factor(subject)54 4.454e-01 2.611e-01 1.706 0.087976 .
#> factor(subject)55 2.685e-01 2.709e-01 0.991 0.321639
#> factor(subject)56 1.125e+00 2.341e-01 4.805 1.55e-06 ***
#> factor(subject)57 2.685e-01 2.709e-01 0.991 0.321639
#> factor(subject)58 -6.006e-01 3.436e-01 -1.748 0.080463 .
#> factor(subject)59 -7.447e-02 2.942e-01 -0.253 0.800210
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> (Dispersion parameter for poisson family taken to be 1)
#>
#> Null deviance: 3185.11 on 117 degrees of freedom
#> Residual deviance: 303.96 on 57 degrees of freedom
#> AIC: 1004.3
#>
#> Number of Fisher Scoring iterations: 5
#>
## IGNORE_RDIFF_END
summary(glmmPQL(y ~ lbase*trt + lage + V4,
random = ~ 1 | subject,
family = poisson, data = epil))
#> iteration 1
#> iteration 2
#> iteration 3
#> iteration 4
#> iteration 5
#> Linear mixed-effects model fit by maximum likelihood
#> Data: epil
#> AIC BIC logLik
#> NA NA NA
#>
#> Random effects:
#> Formula: ~1 | subject
#> (Intercept) Residual
#> StdDev: 0.4449054 1.399283
#>
#> Variance function:
#> Structure: fixed weights
#> Formula: ~invwt
#> Fixed effects: y ~ lbase * trt + lage + V4
#> Value Std.Error DF t-value p-value
#> (Intercept) 1.8706422 0.1056178 176 17.711434 0.0000
#> lbase 0.8828134 0.1293711 54 6.823884 0.0000
#> trtprogabide -0.3103624 0.1491240 54 -2.081236 0.0422
#> lage 0.5375809 0.3465146 54 1.551395 0.1266
#> V4 -0.1610871 0.0773565 176 -2.082398 0.0388
#> lbase:trtprogabide 0.3410064 0.2034656 54 1.675990 0.0995
#> Correlation:
#> (Intr) lbase trtprg lage V4
#> lbase -0.125
#> trtprogabide -0.691 0.088
#> lage -0.103 -0.038 0.088
#> V4 -0.162 0.000 0.000 0.000
#> lbase:trtprogabide 0.055 -0.645 -0.184 0.267 0.000
#>
#> Standardized Within-Group Residuals:
#> Min Q1 Med Q3 Max
#> -2.13615625 -0.63779485 -0.08342573 0.42100427 4.97927910
#>
#> Number of Observations: 236
#> Number of Groups: 59
summary(glmmPQL(y ~ pred, random = ~1 | subject,
family = poisson, data = epil3))
#> iteration 1
#> iteration 2
#> iteration 3
#> iteration 4
#> iteration 5
#> iteration 6
#> iteration 7
#> iteration 8
#> Linear mixed-effects model fit by maximum likelihood
#> Data: epil3
#> AIC BIC logLik
#> NA NA NA
#>
#> Random effects:
#> Formula: ~1 | subject
#> (Intercept) Residual
#> StdDev: 0.7260474 2.169273
#>
#> Variance function:
#> Structure: fixed weights
#> Formula: ~invwt
#> Fixed effects: y ~ pred
#> Value Std.Error DF t-value p-value
#> (Intercept) 3.214199 0.10574213 58 30.396582 0.0000
#> predplacebo-base 0.112880 0.09997031 57 1.129132 0.2636
#> preddrug-placebo -0.107590 0.13493917 57 -0.797319 0.4286
#> Correlation:
#> (Intr) prdpl-
#> predplacebo-base 0.081
#> preddrug-placebo -0.010 -0.700
#>
#> Standardized Within-Group Residuals:
#> Min Q1 Med Q3 Max
#> -2.0412058 -0.4766737 -0.1992579 0.3179334 2.6501959
#>
#> Number of Observations: 118
#> Number of Groups: 59