> library(foreign)
> fev1 <- read.dta("c:/fev1.dta")
> fev1 <- subset(fev1, id != 197)
> attach(fev1)
> loght <- log(ht)
> logbht <- log(baseht)
> library(nlme)
> model1 <- lme(logfev1 ~ age + loght + baseage + logbht,
+ random= ~ age | id)
> summary(model1)
Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
-4549.882 -4499.528 2283.941
Random effects:
Formula: ~age | id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.11048554 (Intr)
age 0.00707838 -0.553
Residual 0.06023788
Fixed effects: logfev1 ~ age + loght + baseage + logbht
Value Std.Error DF t-value p-value
(Intercept) -0.2883233 0.03871675 1692 -7.44699 0.0000
age 0.0235286 0.00139534 1692 16.86231 0.0000
loght 2.2371984 0.04353724 1692 51.38585 0.0000
baseage -0.0165088 0.00745785 296 -2.21362 0.0276
logbht 0.2182148 0.14552087 296 1.49954 0.1348
Correlation:
(Intr) age loght baseag
age 0.023
loght -0.077 -0.875
baseage -0.822 -0.184 0.180
logbht 0.370 0.239 -0.275 -0.815
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-6.45672792 -0.52534885 0.05351814 0.60114614 2.76671671
Number of Observations: 1993
Number of Groups: 299
Linear Mixed Effects Model (Random Intercept and Slope for Log Height)
> model2 <- lme(logfev1 ~ age + loght + baseage + logbht,
+ random= ~ loght | id)
> summary(model2)
Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
-4571.473 -4521.119 2294.737
Random effects:
Formula: ~loght | id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.11534281 (Intr)
loght 0.26170732 -0.615
Residual 0.05944262
Fixed effects: logfev1 ~ age + loght + baseage + logbht
Value Std.Error DF t-value p-value
(Intercept) -0.2846120 0.03901488 1692 -7.29496 0.0000
age 0.0232695 0.00124736 1692 18.65494 0.0000
loght 2.2523371 0.04613191 1692 48.82384 0.0000
baseage -0.0162974 0.00743908 296 -2.19078 0.0292
logbht 0.1807970 0.14548820 296 1.24269 0.2150
Correlation:
(Intr) age loght baseag
age 0.047
loght -0.092 -0.856
baseage -0.823 -0.176 0.148
logbht 0.359 0.247 -0.262 -0.805
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-6.52864459 -0.51673105 0.07986723 0.57825518 2.74313105
Number of Observations: 1993
Number of Groups: 299
Linear Mixed Effects Model (Random Intercept and Slopes for Age and Log Height)
> model3 <- lme(logfev1~ age + loght + baseage + logbht,
+ random= ~ age + loght | id)
> summary(model3)
Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
-4565.899 -4498.761 2294.950
Random effects:
Formula: ~age + loght | id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.115620822 (Intr) age
age 0.003434276 -0.165
loght 0.282546918 -0.509 -0.373
Residual 0.059292700
Fixed effects: logfev1 ~ age + loght + baseage + logbht
Value Std.Error DF t-value p-value
(Intercept) -0.2856675 0.03902885 1692 -7.31939 0.0000
age 0.0234442 0.00127804 1692 18.34393 0.0000
loght 2.2475695 0.04691935 1692 47.90282 0.0000
baseage -0.0161016 0.00742933 296 -2.16731 0.0310
logbht 0.1774347 0.14536410 296 1.22062 0.2232
Correlation:
(Intr) age loght baseag
age 0.035
loght -0.081 -0.861
baseage -0.824 -0.168 0.143
logbht 0.358 0.241 -0.259 -0.804
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-6.54206466 -0.50396522 0.07743504 0.57004654 2.74956818
Number of Observations: 1993
Number of Groups: 299
> library(foreign)
> fat <- read.dta("c:/fat.dta")
> attach(fat)
> time0 <- time*I(time>=0)
> library(nlme)
> model1 <- lme(pbf ~ time + time0, random= ~ time + time0 | id)
> summary(model1)
Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
6082.401 6131.929 -3031.201
Random effects:
Formula: ~time + time0 | id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 6.778082 (Intr) time
time 1.277122 0.292
time0 1.658240 -0.544 -0.827
Residual 3.077864
Fixed effects: pbf ~ time + time0
Value Std.Error DF t-value p-value
(Intercept) 21.361384 0.5645616 885 37.83712 0.0000
time 0.417113 0.1571558 885 2.65414 0.0081
time0 2.047139 0.2279688 885 8.97991 0.0000
Correlation:
(Intr) time
time 0.351
time0 -0.515 -0.872
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.77422856 -0.59003368 -0.03592489 0.59463057 3.37977984
Number of Observations: 1049
Number of Groups: 162
Linear Mixed Effects Model (Hybrid Model with Exponential Serial Correlation)
> model2 <- lme(pbf ~ time + time0, random= ~ 1 | id,
+ corr=corCAR1(, form= ~ time | id))
> summary(model2)
Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
6031.569 6061.285 -3009.784
Random effects:
Formula: ~1 | id
(Intercept) Residual
StdDev: 5.600917 4.159516
Correlation Structure: Continuous AR(1)
Formula: ~time | id
Parameter estimate(s):
Phi
0.4917189
Fixed effects: pbf ~ time + time0
Value Std.Error DF t-value p-value
(Intercept) 21.29795 0.5390841 885 39.50765 0.0000
time 0.24480 0.1428355 885 1.71386 0.0869
time0 2.15342 0.2299134 885 9.36622 0.0000
Correlation:
(Intr) time
time 0.374
time0 -0.454 -0.861
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.4439299 -0.5956301 -0.0118229 0.5562261 3.6872898
Number of Observations: 1049
Number of Groups: 162
> library(foreign)
> cd4 <- read.dta("c:/cd4.dta")
> attach(cd4)
> trt <- group
> trt[group==4] <- 1
> trt[group < 4] <- 0
> week16 <- (week-16)*I(week>=16)
> trt.week <- I(trt==1)*week
> trt.week16 <- I(trt==1)*week16
> library(nlme)
> model1 <- lme(logcd4 ~ week + week16 + trt.week + trt.week16,
+ random= ~ week + week16 | id)
> summary(model1)
Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
11958.97 12037.25 -5967.483
Random effects:
Formula: ~week + week16 | id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.76534028 (Intr) week
week 0.03038383 0.312
week16 0.03521876 -0.458 -0.858
Residual 0.55331915
Fixed effects: logcd4 ~ week + week16 + trt.week + trt.week16
Value Std.Error DF t-value p-value
(Intercept) 2.9414587 0.025620788 3723 114.80750 0e+00
week -0.0073438 0.001986834 3723 -3.69624 2e-04
week16 -0.0120392 0.003174432 3723 -3.79256 2e-04
trt.week 0.0268521 0.003847188 3723 6.97966 0e+00
trt.week16 -0.0277377 0.006198365 3723 -4.47500 0e+00
Correlation:
(Intr) week week16 trt.wk
week -0.193
week16 0.102 -0.867
trt.week 0.001 -0.497 0.438
trt.week16 -0.001 0.434 -0.507 -0.870
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-4.36206287 -0.41946334 0.02941775 0.47734294 3.83987846
Number of Observations: 5036
Number of Groups: 1309
> model2 <- lme(logcd4 ~ age + sex+ week + week16 + trt.week + trt.week16,
+ random= ~ week + week16 | id)
> summary(model2)
Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
11964.37 12055.69 -5968.183
Random effects:
Formula: ~week + week16 | id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.76252587 (Intr) week
week 0.03039088 0.304
week16 0.03523027 -0.452 -0.858
Residual 0.55333186
Fixed effects: logcd4 ~ age + sex + week + week16 + trt.week + trt.week16
Value Std.Error DF t-value p-value
(Intercept) 2.6455449 0.12800208 3723 20.667984 0.0000
age 0.0099963 0.00301564 1306 3.314818 0.0009
sex -0.0926865 0.07537100 1306 -1.229736 0.2190
week -0.0072961 0.00198632 3723 -3.673193 0.0002
week16 -0.0120626 0.00317422 3723 -3.800174 0.0001
trt.week 0.0267911 0.00384335 3723 6.970769 0.0000
trt.week16 -0.0277246 0.00619658 3723 -4.474179 0.0000
Correlation:
(Intr) age sex week week16 trt.wk
age -0.835
sex -0.423 -0.105
week -0.041 0.003 -0.004
week16 0.022 -0.002 0.000 -0.867
trt.week 0.001 -0.003 0.003 -0.497 0.438
trt.week16 0.001 0.000 -0.002 0.434 -0.507 -0.870
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-4.36499603 -0.41614313 0.03007059 0.47652185 3.83302847
Number of Observations: 5036
Number of Groups: 1309
Linear Mixed Effects Model (Random Intercept and Slope): Predicted Means
> predict(model2)
1 1 1 1 1
3.0640938905 3.0703530031 3.0769664295 2.9924438646 2.8915520947
1 2 2 2 2
2.8082760243 3.3171051588 3.5525866637 3.7880681685 3.6089660701
2 2 3 4 4
3.4115878780 3.1995898453 3.5368596183 3.9121735127 3.9706537224
4 4 5 5 5
4.0379217670 3.4399878551 3.4387406771 3.4514432531 3.4641458291
5 5 5 6 6
3.2379650638 3.0117842985 2.7856035333 2.5609162544 2.4769173145
6 6 6 6 7
2.3879769176 2.2486868581 2.1300513958 2.0661710107 3.0424108127
7 8 8 8 8
3.1139582720 2.7894071016 2.7509469489 2.7117748604 2.4573475901
9 10 11 11 11
3.1598179692 2.3714812574 2.8554854037 2.7933384746 2.5964101588
11 12 12 12 13
2.4057901469 3.4912550654 3.9104157599 3.4668677153 0.7076173357
13 13 14 14 14
0.2614050658 0.1697224778 2.6296530337 2.6169714611 2.6084307840
14 15 15 15 16
2.3839488652 3.6205349678 4.1175619237 3.8734648451 3.2002437808
16 16 16 17 17
3.1746592116 3.1445597603 2.9547125217 2.8898549521 2.6268946449
Linear Mixed Effects Model (Random Intercept and Slope): Empirical BLUPs
> re <- random.effects(model2)
> re
(Intercept) week week16
1 0.147099459 8.122817e-03 2.572282e-05
2 0.285957069 9.940208e-03 -1.523399e-02
3 0.381349634 4.620690e-03 -7.962835e-03
4 0.993481712 1.548332e-02 -3.283986e-02
5 0.526535462 8.883961e-03 -1.779782e-02
6 -0.375755865 -4.233150e-03 7.621791e-03
7 0.038904453 1.524586e-02 -1.441896e-02
8 -0.135305455 2.310573e-03 -1.541699e-03
9 0.186636367 2.261412e-03 -3.897092e-03
10 -0.554599810 -6.719906e-03 1.158041e-02
11 -0.119694436 -7.704689e-03 8.001424e-03
12 0.623961218 6.937565e-03 -1.295757e-02
13 -2.263648969 -4.706236e-02 6.222035e-02
14 -0.336248248 5.484486e-03 -1.605586e-03
15 0.646970511 3.929203e-02 -3.297671e-02
16 0.214483043 3.784524e-03 -4.816865e-03
17 -0.108861160 -8.650447e-03 9.775571e-03
18 0.650685247 1.455136e-02 -1.668981e-02
19 -0.158166718 -1.916455e-03 3.302627e-03
20 -0.792272016 -9.599703e-03 1.654317e-02