> library(foreign)
> ds <- read.dta("tlc.dta")
> ds$baseline <- ds$y0
> tlclong <- reshape(ds, idvar="id", varying=c("y0","y1","y4","y6"),
+ v.names="y", timevar="time", time=1:4, direction="long")
> tlclong <- subset(tlclong, time > 1)
> attach(tlclong)
> week <- time
> week[time==2] <- 1
> week[time==3] <- 4
> week[time==4] <- 6
> time <- time - 1
> week.f <- factor(week, c(1,4,6))
> change <- y - baseline
> cbaseline <- baseline - 26.406
> library(nlme)
> model <- gls(y ~ I(week.f==1) + I(week.f==4) + I(week.f==6) +
+ I(week.f==1 & trt=="Succimer") + I(week.f==4 & trt=="Succimer") +
+ I(week.f==6 & trt=="Succimer"), corr=corSymm(, form= ~ time | id),
+ weights = varIdent(form = ~ 1 | week.f))
> summary(model)
Generalized least squares fit by REML
Model: y ~ I(week.f == 1) + I(week.f == 4) + I(week.f == 6) + I(week.f == 1 & trt == "Succimer") + I(week.f == 4 & trt == "Succimer") + I(week.f == 6 & trt == "Succimer")
Data: NULL
AIC BIC logLik
2451.990 2519.544 -1208.995
Correlation Structure: General
Formula: ~time | id
Parameter estimate(s):
Correlation:
1 2 3
2 0.569
3 0.568 0.775
4 0.575 0.581 0.580
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | week.f
Parameter estimates:
0 1 4 6
1.000000 1.330103 1.374827 1.529615
Coefficients:
Value Std.Error t-value p-value
(Intercept) 26.406000 0.4998908 52.82354 0.0000
I(week.f == 1)TRUE -1.644501 0.7824044 -2.10186 0.0362
I(week.f == 4)TRUE -2.231356 0.8073811 -2.76370 0.0060
I(week.f == 6)TRUE -2.642065 0.8864616 -2.98046 0.0031
I(week.f == 1 & trt == "Succimer")TRUE -11.340998 1.0931205 -10.37488 0.0000
I(week.f == 4 & trt == "Succimer")TRUE -8.765288 1.1312570 -7.74827 0.0000
I(week.f == 6 & trt == "Succimer")TRUE -3.119869 1.2507776 -2.49434 0.0130
Correlation:
(Intr) I(.==1 I(.==4 I(.==6 I=1&t="
I(week.f == 1)TRUE -0.155
I(week.f == 4)TRUE -0.136 0.674
I(week.f == 6)TRUE -0.068 0.381 0.380
I(week.f == 1 & trt == "Succimer")TRUE 0.000 -0.699 -0.467 -0.265
I(week.f == 4 & trt == "Succimer")TRUE 0.000 -0.466 -0.701 -0.265 0.667
I(week.f == 6 & trt == "Succimer")TRUE 0.000 -0.263 -0.263 -0.705 0.376
I=4&t="
I(week.f == 1)TRUE
I(week.f == 4)TRUE
I(week.f == 6)TRUE
I(week.f == 1 & trt == "Succimer")TRUE
I(week.f == 4 & trt == "Succimer")TRUE
I(week.f == 6 & trt == "Succimer")TRUE 0.375
Standardized residuals:
Min Q1 Med Q3 Max
-2.1636401 -0.7011814 -0.1426534 0.5374840 5.6570302
Residual standard error: 4.998908
Degrees of freedom: 400 total; 393 residual
> model1 <- gls(y ~ I(week.f==1) + I(week.f==4) + I(week.f==6) +
+ I(week.f==1 & trt=="Succimer") + I(week.f==4 & trt=="Succimer") +
+ I(week.f==6 & trt=="Succimer"), corr=corSymm(, form= ~ time | id),
+ weights = varIdent(form = ~ 1 | week.f), method="ML")
> model2 <- gls(y ~ I(week.f==1) + I(week.f==4) + I(week.f==6),
+ corr=corSymm(, form= ~ time | id),
+ weights = varIdent(form = ~ 1 | week.f), method="ML")
> anova(model1,model2)
Model df AIC BIC logLik Test L.Ratio p-value
model1 1 17 2459.440 2527.295 -1212.720
model2 2 14 2529.588 2585.469 -1250.794 1 vs 2 76.14821 <.0001
Analysis of Response Profiles of Changes in Response from Baseline
> library(foreign)
> ds <- read.dta("tlc.dta")
> ds$baseline <- ds$y0
> tlclong <- reshape(ds, idvar="id", varying=c("y0","y1","y4","y6"),
+ v.names="y", timevar="time", time=1:4, direction="long")
> tlclong <- subset(tlclong, time > 1)
> attach(tlclong)
> week <- time
> week[time==2] <- 1
> week[time==3] <- 4
> week[time==4] <- 6
> time <- time - 1
> week.f <- factor(week, c(1,4,6))
> change <- y - baseline
> cbaseline <- baseline - 26.406
> library(nlme)
> model <- gls(change ~ trt*week.f, corr=corSymm(, form= ~ time | id),
+ weights = varIdent(form = ~ 1 | week.f))
> summary(model)
Generalized least squares fit by REML
Model: change ~ trt * week.f
Data: NULL
AIC BIC logLik
1842.915 1887.118 -909.4573
Correlation Structure: General
Formula: ~time | id
Parameter estimate(s):
Correlation:
1 2
2 0.680
3 0.386 0.385
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | week.f
Parameter estimates:
1 4 6
1.000000 1.029037 1.121984
Coefficients:
Value Std.Error t-value p-value
(Intercept) -1.612 0.7919227 -2.035552 0.0427
trtSuccimer -11.406 1.1199478 -10.184403 0.0000
week.f4 -0.590 0.6426966 -0.918007 0.3594
week.f6 -1.014 0.9343128 -1.085290 0.2787
trtSuccimer:week.f4 2.582 0.9089103 2.840764 0.0048
trtSuccimer:week.f6 8.254 1.3213178 6.246794 0.0000
Correlation:
(Intr) trtScc wek.f4 wek.f6 trS:.4
trtSuccimer -0.707
week.f4 -0.369 0.261
week.f6 -0.480 0.340 0.325
trtSuccimer:week.f4 0.261 -0.369 -0.707 -0.230
trtSuccimer:week.f6 0.340 -0.480 -0.230 -0.707 0.325
Standardized residuals:
Min Q1 Med Q3 Max
-3.32746749 -0.45341726 -0.01492447 0.49181724 5.72640936
Residual standard error: 5.599739
Degrees of freedom: 300 total; 294 residual
> model1 <- gls(change ~ trt*week.f, corr=corSymm(, form= ~ time | id),
+ weights = varIdent(form = ~ 1 | week.f),method="ML")
> model2 <- gls(change ~ week.f, corr=corSymm(, form= ~ time | id),
+ weights = varIdent(form = ~ 1 | week.f),method="ML")
> anova(model1,model2)
Model df AIC BIC logLik Test L.Ratio p-value
model1 1 12 1850.023 1894.469 -913.0117
model2 2 9 1918.211 1951.545 -950.1055 1 vs 2 74.18778 <.0001
Analysis of Response Profiles of Adjusted Changes in Response from Baseline
> model <- gls(change ~ cbaseline + trt*week.f,
+ corr=corSymm(, form= ~ time | id),
+ weights = varIdent(form = ~ 1 | week.f))
> summary(model)
Generalized least squares fit by REML
Model: change ~ cbaseline + trt * week.f
Data: NULL
AIC BIC logLik
1843.566 1891.408 -908.7828
Correlation Structure: General
Formula: ~time | id
Parameter estimate(s):
Correlation:
1 2
2 0.669
3 0.377 0.377
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | week.f
Parameter estimates:
1 4 6
1.000000 1.033903 1.144265
Coefficients:
Value Std.Error t-value p-value
(Intercept) -1.638195 0.7766433 -2.109327 0.0358
cbaseline -0.195482 0.0938963 -2.081889 0.0382
trtSuccimer -11.353611 1.0984835 -10.335713 0.0000
week.f4 -0.590000 0.6427037 -0.917997 0.3594
week.f6 -1.014000 0.9343107 -1.085292 0.2787
trtSuccimer:week.f4 2.582000 0.9089202 2.840733 0.0048
trtSuccimer:week.f6 8.254000 1.3213149 6.246808 0.0000
Correlation:
(Intr) cbasln trtScc wek.f4 wek.f6 trS:.4
cbaseline 0.016
trtSuccimer -0.707 -0.023
week.f4 -0.372 0.000 0.263
week.f6 -0.473 0.000 0.334 0.325
trtSuccimer:week.f4 0.263 0.000 -0.372 -0.707 -0.230
trtSuccimer:week.f6 0.334 0.000 -0.473 -0.230 -0.707 0.325
Standardized residuals:
Min Q1 Med Q3 Max
-3.01806838 -0.50986122 -0.06341362 0.41107468 5.94888816
Residual standard error: 5.490976
Degrees of freedom: 300 total; 293 residual
> model1 <- gls(change ~ cbaseline + trt*week.f,
+ corr=corSymm(, form= ~ time | id),
+ weights = varIdent(form = ~ 1 | week.f),method="ML")
> model2 <- gls(change ~ cbaseline + week.f,
+ corr=corSymm(, form= ~ time | id),
+ weights = varIdent(form = ~ 1 | week.f),method="ML")
> anova(model1,model2)
Model df AIC BIC logLik Test L.Ratio p-value
model1 1 13 1847.689 1895.839 -910.8447
model2 2 10 1917.839 1954.877 -948.9195 1 vs 2 76.14967 <.0001