> library(foreign)
> ds <- read.dta("smoking.dta")
> attach(ds)
> year <- time
> year[time==0] <- 1
> year[time==3] <- 2
> year[time==6] <- 3
> year[time==9] <- 4
> year[time==12] <- 5
> year[time==15] <- 6
> year[time==19] <- 7
> library(nlme)
> model <- gls(fev1 ~ smoker*time, corr=corSymm(, form= ~ year | id),
+ weights = varIdent(form = ~ 1 | year))
> summary(model)
Generalized least squares fit by REML
Model: fev1 ~ smoker * time
Data: NULL
AIC BIC logLik
330.0506 478.6102 -133.0253
Correlation Structure: General
Formula: ~year | id
Parameter estimate(s):
Correlation:
1 2 3 4 5 6
2 0.863
3 0.846 0.888
4 0.838 0.833 0.833
5 0.855 0.862 0.890 0.886
6 0.839 0.874 0.868 0.876 0.932
7 0.831 0.824 0.834 0.840 0.858 0.893
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | year
Parameter estimates:
1 2 3 4 6 7 5
1.0000000 0.9779830 0.9509433 0.9550638 0.9679868 0.9198133 0.9387381
Coefficients:
Value Std.Error t-value p-value
(Intercept) 3.507313 0.10037735 34.94128 0.0000
smoker -0.261700 0.11509873 -2.27370 0.0233
time -0.033224 0.00306635 -10.83512 0.0000
smoker:time -0.004998 0.00352539 -1.41784 0.1566
Correlation:
(Intr) smoker time
smoker -0.872
time -0.369 0.322
smoker:time 0.321 -0.364 -0.870
Standardized residuals:
Min Q1 Med Q3 Max
-2.962521431 -0.676944345 -0.002825979 0.675168413 3.202172472
Residual standard error: 0.595934
Degrees of freedom: 771 total; 767 residual
Linear Trend Model (ML Estimation)
> model1 <- gls(fev1 ~ smoker*time, corr=corSymm(, form= ~ year | id),
+ weights = varIdent(form = ~ 1 | year),method="ML" )
> summary(model1)
Generalized least squares fit by maximum likelihood
Model: fev1 ~ smoker * time
Data: NULL
AIC BIC logLik
302.4609 451.1869 -119.2305
Correlation Structure: General
Formula: ~year | id
Parameter estimate(s):
Correlation:
1 2 3 4 5 6
2 0.861
3 0.844 0.886
4 0.836 0.831 0.831
5 0.854 0.861 0.889 0.884
6 0.839 0.873 0.867 0.875 0.931
7 0.831 0.823 0.832 0.838 0.856 0.891
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | year
Parameter estimates:
1 2 3 4 6 7 5
1.0000000 0.9779636 0.9508766 0.9551056 0.9682082 0.9189570 0.9386682
Coefficients:
Value Std.Error t-value p-value
(Intercept) 3.507419 0.09988395 35.11494 0.0000
smoker -0.261739 0.11453237 -2.28529 0.0226
time -0.033226 0.00304990 -10.89424 0.0000
smoker:time -0.005005 0.00350642 -1.42746 0.1539
Correlation:
(Intr) smoker time
smoker -0.872
time -0.369 0.322
smoker:time 0.321 -0.364 -0.870
Standardized residuals:
Min Q1 Med Q3 Max
-2.981881123 -0.681580993 -0.002822034 0.679638296 3.223786402
Residual standard error: 0.5919557
Degrees of freedom: 771 total; 767 residual
Quadratic Trend Model (ML Estimation)
> model2 <- gls(fev1 ~ smoker*time + smoker*timesqr,
+ corr=corSymm(, form= ~ year | id),
+ weights = varIdent(form = ~ 1 | year),method="ML")
> summary(model2)
Generalized least squares fit by maximum likelihood
Model: fev1 ~ smoker * time + smoker * timesqr
Data: NULL
AIC BIC logLik
305.1824 463.2038 -118.5912
Correlation Structure: General
Formula: ~year | id
Parameter estimate(s):
Correlation:
1 2 3 4 5 6
2 0.860
3 0.848 0.888
4 0.834 0.830 0.830
5 0.855 0.861 0.888 0.884
6 0.836 0.873 0.867 0.874 0.932
7 0.831 0.822 0.835 0.836 0.857 0.891
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | year
Parameter estimates:
1 2 3 4 6 7 5
1.0000000 0.9778155 0.9501207 0.9556550 0.9691949 0.9189935 0.9391731
Coefficients:
Value Std.Error t-value p-value
(Intercept) 3.531198 0.10380634 34.01717 0.0000
smoker -0.297785 0.11883818 -2.50581 0.0124
time -0.041339 0.01008819 -4.09780 0.0000
timesqr 0.000411 0.00048706 0.84305 0.3995
smoker:time 0.007550 0.01152722 0.65498 0.5127
smoker:timesqr -0.000642 0.00056132 -1.14413 0.2529
Correlation:
(Intr) smoker time timsqr smokr:tm
smoker -0.874
time -0.366 0.319
timesqr 0.271 -0.237 -0.953
smoker:time 0.320 -0.359 -0.875 0.834
smoker:timesqr -0.235 0.265 0.827 -0.868 -0.952
Standardized residuals:
Min Q1 Med Q3 Max
-2.97161572 -0.68806223 0.01034355 0.66356295 3.25067668
Residual standard error: 0.5914762
Degrees of freedom: 771 total; 765 residual
> anova(model2,model1)
Model df AIC BIC logLik Test L.Ratio p-value
model2 1 34 305.1824 463.2038 -118.5912
model1 2 32 302.4609 451.1869 -119.2305 1 vs 2 1.278513 0.5277
> library(foreign)
> ds <- read.dta("c:/tlc.dta")
> tlclong <- reshape(ds, idvar="id", varying=c("y0","y1","y4","y6"),
+ v.names="y", timevar="time", time=1:4, direction="long")
> attach(tlclong)
> week <- time
> week[time==1] <- 0
> week[time==2] <- 1
> week[time==3] <- 4
> week[time==4] <- 6
> week1 <- (week-1)*I(week>=1)
> trt.week <- week*I(trt=="Succimer")
> trt.week1 <- week1*I(trt=="Succimer")
> library(nlme)
> model <- gls(y ~ week + week1 + trt.week + trt.week1,
+ corr=corSymm(, form= ~ time | id),
+ weights = varIdent(form = ~ 1 | time))
> summary(model)
Generalized least squares fit by REML
Model: y ~ week + week1 + trt.week + trt.week1
Data: NULL
AIC BIC logLik
2467.452 2527.136 -1218.726
Correlation Structure: General
Formula: ~time | id
Parameter estimate(s):
Correlation:
1 2 3
2 0.569
3 0.560 0.768
4 0.574 0.576 0.553
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | time
Parameter estimates:
1 2 3 4
1.000000 1.329903 1.385011 1.544057
Coefficients:
Value Std.Error t-value p-value
(Intercept) 26.342207 0.4991198 52.77733 0.0000
week -1.629603 0.7817592 -2.08453 0.0378
week1 1.430494 0.8777585 1.62971 0.1040
trt.week -11.249985 1.0924522 -10.29792 0.0000
trt.week1 12.582249 1.2278544 10.24735 0.0000
Correlation:
(Intr) week week1 trt.wk
week -0.154
week1 0.147 -0.988
trt.week 0.000 -0.699 0.691
trt.week1 0.000 0.690 -0.699 -0.987
Standardized residuals:
Min Q1 Med Q3 Max
-2.0020271 -0.6888161 -0.1136309 0.5520751 5.8000806
Residual standard error: 4.999256
Degrees of freedom: 400 total; 395 residual