> library(foreign)
> ds <- read.dta("c:/exercise.dta")
> exlong <- reshape(ds, idvar="id",
+ varying=c("y0","y2","y4","y6","y8","y10","y12"), v.names="y",
+ timevar="time", time=0:6, direction="long")
> exlong <- subset(exlong, time!=1 & time!=5)
> attach(exlong)
> day <- time*2
> day.f <- factor(day, c(0,4,6,8,12))
> group.f <- factor(group, c(1,2))
> newtime <- time
> newtime[time==0] <- 1
> newtime[time==2] <- 2
> newtime[time==3] <- 3
> newtime[time==4] <- 4
> newtime[time==6] <- 5
> library(nlme)
> model1 <- gls(y ~ group.f*day.f, na.action=na.omit,
+ corr=corSymm(, form= ~ newtime | id),
+ weights = varIdent(form = ~ 1 | newtime))
> summary(model1)
Generalized least squares fit by REML
Model: y ~ group.f * day.f
Data: NULL
AIC BIC logLik
647.3399 724.6837 -298.6699
Correlation Structure: General
Formula: ~newtime | id
Parameter estimate(s):
Correlation:
1 2 3 4
2 0.924
3 0.885 0.960
4 0.844 0.949 0.958
5 0.810 0.902 0.911 0.939
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | newtime
Parameter estimates:
1 2 3 4 5
1.000000 1.139327 1.049132 1.202937 1.200952
Coefficients:
Value Std.Error t-value p-value
(Intercept) 79.68750 0.7773470 102.51214 0.0000
group.f2 1.36012 1.0318252 1.31817 0.1893
day.f4 1.12500 0.3416830 3.29253 0.0012
day.f6 1.36017 0.3873739 3.51126 0.0006
day.f8 1.58458 0.5048462 3.13873 0.0020
day.f12 1.62356 0.5537518 2.93193 0.0039
group.f2:day.f4 -0.16916 0.4548986 -0.37186 0.7105
group.f2:day.f6 0.21126 0.5123690 0.41231 0.6807
group.f2:day.f8 -0.13092 0.6714472 -0.19498 0.8457
group.f2:day.f12 0.32052 0.7532369 0.42552 0.6710
Correlation:
(Intr) grp.f2 day.f4 day.f6 day.f8 dy.f12 g.2:.4 g.2:.6 g.2:.8
group.f2 -0.753
day.f4 0.119 -0.090
day.f6 -0.144 0.109 0.760
day.f8 0.023 -0.017 0.819 0.820
day.f12 -0.038 0.029 0.665 0.696 0.798
group.f2:day.f4 -0.090 0.119 -0.751 -0.571 -0.615 -0.500
group.f2:day.f6 0.109 -0.145 -0.575 -0.756 -0.620 -0.526 0.761
group.f2:day.f8 -0.017 0.023 -0.616 -0.617 -0.752 -0.600 0.815 0.822
group.f2:day.f12 0.028 -0.037 -0.489 -0.511 -0.586 -0.735 0.647 0.681 0.779
Standardized residuals:
Min Q1 Med Q3 Max
-2.26656151 -0.65852801 -0.01531460 0.59847743 2.53952400
Residual standard error: 3.109388
Degrees of freedom: 173 total; 163 residual
Autoregressive Covariance (REML Estimation)
> model2 <- gls(y ~ group.f*day.f, na.action=na.omit,
+ corr=corAR1(, form= ~ newtime | id))
> summary(model2)
Generalized least squares fit by REML
Model: y ~ group.f * day.f
Data: NULL
AIC BIC logLik
645.0718 682.1968 -310.5359
Correlation Structure: ARMA(1,0)
Formula: ~newtime | id
Parameter estimate(s):
Phi1
0.9401763
Coefficients:
Value Std.Error t-value p-value
(Intercept) 79.68750 0.8612227 92.52833 0.0000
group.f2 1.36012 1.1431591 1.18979 0.2359
day.f4 1.12500 0.2978976 3.77647 0.0002
day.f6 1.31197 0.4186010 3.13418 0.0020
day.f8 1.53493 0.5036625 3.04754 0.0027
day.f12 1.61594 0.5745615 2.81247 0.0055
group.f2:day.f4 -0.21440 0.3976147 -0.53921 0.5905
group.f2:day.f6 0.25946 0.5535425 0.46872 0.6399
group.f2:day.f8 -0.04774 0.6707040 -0.07118 0.9433
group.f2:day.f12 0.35890 0.7788361 0.46082 0.6455
Correlation:
(Intr) grp.f2 day.f4 day.f6 day.f8 dy.f12 g.2:.4 g.2:.6 g.2:.8
group.f2 -0.753
day.f4 -0.173 0.130
day.f6 -0.239 0.180 0.690
day.f8 -0.289 0.218 0.557 0.792
day.f12 -0.328 0.247 0.475 0.674 0.839
group.f2:day.f4 0.130 -0.172 -0.749 -0.517 -0.417 -0.356
group.f2:day.f6 0.181 -0.240 -0.522 -0.756 -0.599 -0.510 0.689
group.f2:day.f8 0.217 -0.288 -0.418 -0.595 -0.751 -0.630 0.552 0.792
group.f2:day.f12 0.242 -0.321 -0.350 -0.497 -0.619 -0.738 0.462 0.663 0.824
Standardized residuals:
Min Q1 Med Q3 Max
-2.18724319 -0.59439302 -0.01382309 0.58072348 2.62469204
Residual standard error: 3.444891
Degrees of freedom: 173 total; 163 residual
Exponential Covariance (REML Estimation)
> model3 <- gls(y ~ group.f*day.f, na.action=na.omit,
+ corr=corExp(, form= ~ day | id))
> summary(model3)
Generalized least squares fit by REML
Model: y ~ group.f * day.f
Data: NULL
AIC BIC logLik
642.5459 679.6709 -309.2729
Correlation Structure: Exponential spatial correlation
Formula: ~day | id
Parameter estimate(s):
range
46.13188
Coefficients:
Value Std.Error t-value p-value
(Intercept) 79.68750 0.8614996 92.49859 0.0000
group.f2 1.36012 1.1435267 1.18941 0.2360
day.f4 1.12500 0.3511182 3.20405 0.0016
day.f6 1.31042 0.4279933 3.06178 0.0026
day.f8 1.57188 0.4890739 3.21399 0.0016
day.f12 1.60783 0.5898073 2.72603 0.0071
group.f2:day.f4 -0.21083 0.4678108 -0.45067 0.6528
group.f2:day.f6 0.26101 0.5666630 0.46061 0.6457
group.f2:day.f8 -0.07356 0.6507037 -0.11305 0.9101
group.f2:day.f12 0.34980 0.8031588 0.43554 0.6637
Correlation:
(Intr) grp.f2 day.f4 day.f6 day.f8 dy.f12 g.2:.4 g.2:.6 g.2:.8
group.f2 -0.753
day.f4 -0.204 0.154
day.f6 -0.245 0.185 0.803
day.f8 -0.280 0.211 0.688 0.846
day.f12 -0.335 0.252 0.548 0.673 0.785
group.f2:day.f4 0.153 -0.203 -0.751 -0.603 -0.516 -0.411
group.f2:day.f6 0.185 -0.246 -0.606 -0.755 -0.639 -0.509 0.802
group.f2:day.f8 0.211 -0.280 -0.517 -0.636 -0.752 -0.590 0.684 0.846
group.f2:day.f12 0.246 -0.326 -0.402 -0.495 -0.576 -0.734 0.532 0.658 0.766
Standardized residuals:
Min Q1 Med Q3 Max
-2.18976671 -0.59420196 -0.01381865 0.58098711 2.62281222
Residual standard error: 3.445998
Degrees of freedom: 173 total; 163 residual
> anova(model1,model2)
Model df AIC BIC logLik Test L.Ratio p-value
model1 1 25 647.3399 724.6837 -298.6699
model2 2 12 645.0718 682.1968 -310.5359 1 vs 2 23.73189 0.0337
> anova(model1,model3)
Model df AIC BIC logLik Test L.Ratio p-value
model1 1 25 647.3399 724.6837 -298.6699
model3 2 12 642.5459 679.6709 -309.2729 1 vs 2 21.20597 0.069
> anova(model2,model3)
Model df AIC BIC logLik
model2 1 12 645.0718 682.1968 -310.5359
model3 2 12 642.5459 679.6709 -309.2729