8452 HW1
Aparajita Sur
10/12/2020
Problem 2 Simulations
library(mvtnorm)
library(gee)
library(lattice)
#Initial parameters
set.seed(1996)
n = 10 #observations per subject
m = 25 #number of subjects
beta = c(10,0) #Beta vector
t = 1:10 #times tij
Mean.model = beta[1] + beta[2] * t
### Functions to generate data
##Random Intercepts
RI.data = function(sigma, tau, n=10, t.ij){
b0i = rnorm(m, mean = 0, sd = tau)
b0i = rep(b0i, rep(n,m))
epsilon = rnorm((n*m), mean = 0, sd = sigma )
Y = rep(Mean.model, m) + b0i + epsilon
return( data.frame(id=rep(1:m,each=n),
y=Y,
t.ij=rep(t.ij,m)))}
##Random Intercepts and Slopes
RIS.data = function(sigma,D, n=10, t.ij){
bi = rmvnorm(m,mean = c(0, 0),sigma = D )
b0i = rep(bi[,1], rep(n, m))
b1i = rep(bi[,2], rep(n, m))
epsilon = rnorm((n*m), mean = 0, sd = sigma )
Y = rep( Mean.model, m ) + b0i + ( b1i * rep(t, m) ) + epsilon
return( data.frame(id=rep(1:m,each=n),
y=Y,
t.ij=rep(t.ij,m)))}
## Serial Correlation
SC.data = function( sigma, tau, rho, n=10, t.ij){
Sigma = matrix( tau^2, nrow = n, ncol = n ) + diag( rep(sigma^2, n) )
for( i in 1:n ){
for( j in 1:n ){
Sigma[i,j] = Sigma[i,j] * rho^(abs(i-j))}}
return(data.frame( id=rep(1:m,each=n),y=as.vector(t(rmvnorm( m, Mean.model,
sigma = Sigma))),t.ij=rep(t.ij,m)))}#Random intercept model plots
#par(mfrow=c(1,3))
data.ri= RI.data(sigma=1,tau=1,n=10,t.ij=t)
ri.11=plot(y~t.ij,data=data.ri,xlab="tij",ylab="Yi",col="gray50",
main="Sigma=1, Tau=1")
for (i in 1:25){lines(data.ri$t.ij[data.ri$id==i],
data.ri$y[data.ri$id==i],lwd=1.5)}
with (data.ri, {lines(supsmu(t.ij,y),lty=5,col="red",lwd=2)})
data.ri= RI.data(sigma=1,tau=2,n=10,t.ij=t)
ri.12=plot(y~t.ij,data=data.ri,xlab="tij",ylab="Yi",col="gray50",
main="Sigma=1, Tau=2")
for (i in 1:25){lines(data.ri$t.ij[data.ri$id==i],
data.ri$y[data.ri$id==i],lwd=1.5)}
with (data.ri, {lines(supsmu(t.ij,y),lty=5,col="red",lwd=2)})
data.ri= RI.data(sigma=1,tau=5,n=10,t.ij=t)
ri.15=plot(y~t.ij,data=data.ri,xlab="tij",ylab="Yi",col="gray50",
main= "Sigma=1, Tau=5")
for (i in 1:25){lines(data.ri$t.ij[data.ri$id==i],
data.ri$y[data.ri$id==i],lwd=1.5)}
with (data.ri, {lines(supsmu(t.ij,y),lty=5,col="red",lwd=2)})
\(Var(Y_{i})\) for the random intercept model is \(Var(b_{0,1}) + Var(\epsilon_{ij})=\tau^2 + \sigma^2\). \(Cov(Y_{ij},Y_{ik}) = Cov(b_{0,i}+\epsilon_{ij},b_{0,i}+\epsilon_{ij})=Cov(b_{0,i},b_{0,i})=\tau^2\). So \(Cov(Y_{i})\) has \(\tau^2 + \sigma^2\) on the diagonals and \(\tau^2\) on the off diagonals.
#Random intercept and slope model plots
#par(mfrow=c(1,3))
data.ris= RIS.data(sigma=1,D=matrix(c(2,0,0,2),nrow=2,ncol=2),n=10,t.ij=t)
ris.1=plot(y~t.ij,data=data.ri,xlab="tij",ylab="Yi",col="gray50",
main="RIS 1")
for (i in 1:25){lines(data.ris$t.ij[data.ris$id==i],data.ris$y[data.ris$id==i],lwd=1.5)}
with (data.ris, {lines(supsmu(t.ij,y),lty=5,col="red",lwd=2)})
data.ris= RIS.data(sigma=1,D=matrix(c(2,-2,-2,2),nrow=2,ncol=2),n=10,t.ij=t)
ris.2=plot(y~t.ij,data=data.ri,xlab="tij",ylab="Yi",col="gray50",
main="RIS 2")
for (i in 1:25){lines(data.ris$t.ij[data.ris$id==i],data.ris$y[data.ris$id==i],lwd=1.5)}
with (data.ris, {lines(supsmu(t.ij,y),lty=5,col="red",lwd=2)})
data.ris= RIS.data(sigma=1,D=matrix(c(.5,0,0,.4),nrow=2,ncol=2),n=10,t.ij=t)
ris.3=plot(y~t.ij,data=data.ri,xlab="tij",ylab="Yi",col="gray50",
main="RIS 3")
for (i in 1:25){lines(data.ris$t.ij[data.ris$id==i],data.ris$y[data.ris$id==i],lwd=1.5)}
with (data.ris, {lines(supsmu(t.ij,y),lty=5,col="red",lwd=2)})
\(Var(Y_{i})= Var(b_{0,1})+2t_{ij}Cov(b_{0i},b_{1i})+t_{ij}^2Var(b_{1i})+Var(\epsilon_{ij}) = Var(b_{0,1})+2t_{ij}Cov(b_{0i},b_{1i})+t_{ij}^2Var(b_{1i})+\sigma^2\), where \(Var(b_{0,1})\) and \(Var(b_{1i})\) are the diagonals of D (\(d_{11}, d_{22}\)) and the covariance are the off-diagonals (\(d_{12}\)). The \(Cov(Y_{ij},Y_{ik})=d_{11}+(t_{ij}+t_{ik})d_{12}+t_{ij}t_{ik}d_{22}\).
####### Serial Correlation
data.sc= SC.data(sigma=1, tau=2, rho=.7, n=10, t.ij=t)
SC.12.7=plot(y~t.ij,data=data.sc,xlab="tij",ylab="Yi",col="gray50",
main="s=1,tau=2,rho=.7")
for (i in 1:25){lines(data.sc$t.ij[data.sc$id==i],data.sc$y[data.sc$id==i],lwd=1.5)}
with (data.sc, {lines(supsmu(t.ij,y),lty=5,col="red",lwd=2)})
data.sc= SC.data(sigma=1, tau=2, rho=.9, n=10, t.ij=t)
SC.12.9=plot(y~t.ij,data=data.sc,xlab="tij",ylab="Yi",col="gray50",
main="s=1,tau=2,rho=.7")
for (i in 1:25){lines(data.sc$t.ij[data.sc$id==i],data.sc$y[data.sc$id==i],lwd=1.5)}
with (data.sc, {lines(supsmu(t.ij,y),lty=5,col="red",lwd=2)})
data.sc= SC.data(sigma=2, tau=2, rho=.9, n=10, t.ij=t)
SC.22.9=plot(y~t.ij,data=data.sc,xlab="tij",ylab="Yi",col="gray50",
main="s=1,tau=2,rho=.7")
for (i in 1:25){lines(data.sc$t.ij[data.sc$id==i],data.sc$y[data.sc$id==i],lwd=1.5)}
with (data.sc, {lines(supsmu(t.ij,y),lty=5,col="red",lwd=2)})
##3. FEV1 data analysis
#Open libraries and data
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(mvtnorm)
library(gee)
library(lattice)
library(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse#Open and wrangle data
fev1 = read.table("~/Biostat 3rd Year/Fall - Longitudinal/8452 Datasets/fev1.dat",
quote="\"", comment.char="")
names(fev1) = c("id","height","age","base.height","base.age","l.fev1")
fev1$t.ij <- fev1$age - fev1$base.age #time
fev1$id=factor(fev1$id)
fev1$c.age = fev1$age-fev1$base.age #change in age variable Basic scatter plots and nonparametric curve to highlight response over time
summary(fev1)## id height age base.height
## 18 : 12 Min. :1.110 Min. : 6.434 Min. :1.110
## 23 : 12 1st Qu.:1.370 1st Qu.: 9.717 1st Qu.:1.220
## 35 : 12 Median :1.540 Median :12.595 Median :1.260
## 38 : 12 Mean :1.497 Mean :12.566 Mean :1.276
## 39 : 12 3rd Qu.:1.620 3rd Qu.:15.366 3rd Qu.:1.320
## 40 : 12 Max. :1.790 Max. :18.691 Max. :1.720
## (Other):1922
## base.age l.fev1 t.ij c.age
## Min. : 6.434 Min. :-0.6932 Min. : 0.000 Min. : 0.000
## 1st Qu.: 7.136 1st Qu.: 0.5481 1st Qu.: 1.074 1st Qu.: 1.074
## Median : 7.781 Median : 0.8671 Median : 4.108 Median : 4.108
## Mean : 8.030 Mean : 0.8152 Mean : 4.536 Mean : 4.536
## 3rd Qu.: 8.449 3rd Qu.: 1.0978 3rd Qu.: 7.140 3rd Qu.: 7.140
## Max. :14.067 Max. : 1.5953 Max. :11.061 Max. :11.061
## #summary(fev1$base.age)
hist(fev1$base.age)
hist(fev1$l.fev1)
summary(fev1$age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.434 9.717 12.595 12.566 15.366 18.691#######Basic scatter plots##########
par(mfrow=c(1,3))
plot(l.fev1~base.age, data=fev1,pch=".",xlab="Baseline Age",
ylab="Log(FEV1)")
with (fev1,{lines(supsmu(base.age,l.fev1),lty=5,col="red",lwd=2)})
#Among 7-9 year olds
subset=which(fev1$base.age<10)
plot(fev1$l.fev1[subset]~fev1$base.age[subset],pch=".",xlab="Baseline Age under 10",
ylab="Log(FEV1)")
with (fev1,{lines(supsmu(base.age[subset],l.fev1[subset]),lty=5,col="red",lwd=2)})
plot(l.fev1~c.age, data=fev1,pch=".",xlab="Change in Age",
ylab="Log(FEV1)")
with (fev1,{lines(supsmu(c.age,l.fev1),lty=5,col="red",lwd=2)})
plot(l.fev1[subset]~c.age[subset], data=fev1,pch=".",xlab="Change in Age for girls starting at age 10",
ylab="Log(FEV1)")
with (fev1,{lines(supsmu(c.age[subset],l.fev1[subset]),lty=5,col="red",lwd=2)})
plot(l.fev1~height, data=fev1,pch=".",xlab="Height",
ylab="Log(FEV1)")
plot(l.fev1~t.ij, data=fev1,pch=".",xlab="Annual examination (tij)",
ylab="Log(FEV1)")
with (fev1,{lines(supsmu(t.ij,l.fev1),lty=5,col="red",lwd=2)})
plot(height~age,data=fev1)
Spaghetti Plots (all indv, age groups and height groups)
###### Spaghetti plot #####
xyplot (l.fev1 ~ t.ij, data = fev1, type = "b", group = id,
xlab = "Annual examiniation (tij)",
ylab = "Log(FEV1)",
col.line = "gray20")
#Plot by baseline age
xyplot (l.fev1 ~ base.age | equal.count (base.age, 4), data = fev1,
type = "l", group = id,
xlab = "Annual examination (tij)",
col.line = "gray20",
ylab = "Log(FEV1)",
strip = strip.custom (var.name = "Baseline Age"))
#Plot by age
xyplot (l.fev1 ~ age | equal.count (age, 4), data = fev1,
type = "l", group = id,
xlab = "Annual examination (tij)",
col.line = "gray20",
ylab = "Log(FEV1)",
strip = strip.custom (var.name = "Age"))
#Plot by height
xyplot (l.fev1 ~ height | equal.count (height, 4), data = fev1,
type = "l", group = id,
xlab = "Annual examination (tij)",
col.line = "gray20",
ylab = "Log(FEV1)",
strip = strip.custom (var.name = "Height"))
The maximum followup period was 11 years, while the median followup was near 4.The baseline age ranged from 6 to 14 years old while the majority of the subjects were between 7-8 years old at study entry. At baseline, the age did not have an effect on log(fev1) but as the change in age increased the log(fev1) slightly increased, the same trend as follow-up time. This could be due to the small variability among baseline ages, but from this data it seems that there is a stronger longitudinal effect compared to the cross-sectional effect of age. Height seems correlated with log(fev1) as well, so a model including both age and height might be optimal.
Exploring the correlation structure
#Source: this code is from Professor Luo's lecture notes
fev1.lm <- lm (l.fev1 ~ t.ij, data = fev1)
fev1$lmres <- resid (fev1.lm)
fev1$roundyr <- round (fev1$t.ij)
## Reshape the data to wide format
fev1w <- reshape (fev1[,c("id", "lmres", "roundyr")],
direction = "wide",
v.names = "lmres", timevar = "roundyr",
idvar = "id")
#Put histograms on the diagonal
panel.hist <- function(x, ...) {
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks;
nB <- length(breaks)
y <- h$counts;
y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col="cyan", ...)}
## Put (absolute) correlations on the upper panel,
#w/ size prop. to correlation.
panel.cor <- function(x, y, digits=2, prefix="", cex.cor)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs (cor(x, y, use = "pairwise.complete.obs"))
txt <- format(c(r, 0.123456789), digits=digits)[1]
txt <- paste(prefix, txt, sep="")
if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex * r)
}
pairs (fev1w[,c(5,2,3,6:9)], upper.panel = panel.cor,
diag.panel = panel.hist)
From the correlation histogram, the correlations for farther observations are only slight smaller than closer observations and the correlations are constant along the diagonal. There may be evidence of strong within subject correlation consistently over time which would make a random effects model appealing.
Variogram
#Source: Professor Luo's lecture code
lda.vg <- function (id, res, time, plot = TRUE, ...) {
vv <- tapply (res, id,
function (x) outer(x, x, function(x, y) (x - y)^2/2))
v <- unlist (lapply (vv, function (x) x[lower.tri (x)]))
uu <- tapply (time, id,
function (x) outer (x, x, function (x, y) (x - y)))
u <- unlist (lapply (uu, function (x) x[lower.tri (x)]))
if (plot) {
vg.loess <- loess.smooth (u, v, family = "gaussian")
plot (v ~ u, pch = ".", col = "gray50", ...)
lines (vg.loess, lty = 1)
abline (h = var (res), lty = 2)
}
invisible (data.frame (v = v, u = u))
}
lda.vg (fev1$id, fev1$lmres, fev1$t.ij,
ylim = c(0, .05), ylab = "", xlab = "",main="Variogram of Log(FEV1) residuals")
There is a smooth rise which levels off at the 6 year followup after which it decreases. This trend may be due to the lack of data for longer follow-up times. There does not seem to be a large autocorrelation effect.
b. Modeling Log(FEV1) with age and change in age as covariates
#OLS model
fev1.age.lm = lm(l.fev1~base.age+c.age,data=fev1)
summary(fev1.age.lm)##
## Call:
## lm(formula = l.fev1 ~ base.age + c.age, data = fev1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.21787 -0.10089 0.00294 0.11044 0.58188
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.312591 0.026447 -11.82 <2e-16 ***
## base.age 0.091648 0.003077 29.79 <2e-16 ***
## c.age 0.086395 0.001106 78.13 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.162 on 1991 degrees of freedom
## Multiple R-squared: 0.7599, Adjusted R-squared: 0.7597
## F-statistic: 3151 on 2 and 1991 DF, p-value: < 2.2e-16#GLS with baseline age and change in age as covariates
fev1.ind=gls(l.fev1~base.age+c.age,data=fev1,method="ML")
summary(fev1.ind) #same as OLS as expected## Generalized least squares fit by maximum likelihood
## Model: l.fev1 ~ base.age + c.age
## Data: fev1
## AIC BIC logLik
## -1594.15 -1571.759 801.0751
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.31259141 0.026446775 -11.81964 0
## base.age 0.09164756 0.003076561 29.78896 0
## c.age 0.08639534 0.001105816 78.12811 0
##
## Correlation:
## (Intr) base.g
## base.age -0.973
## c.age -0.381 0.205
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -7.52163837 -0.62308337 0.01815619 0.68207352 3.59375231
##
## Residual standard error: 0.1619154
## Degrees of freedom: 1994 total; 1991 residualanova(fev1.ind)## Denom. DF: 1991
## numDF F-value p-value
## (Intercept) 1 50474.65 <.0001
## base.age 1 198.54 <.0001
## c.age 1 6104.00 <.0001fev1.ar1=gls(l.fev1~base.age+c.age,data=fev1,correlation = corAR1 (form = ~ 1 | id),method="ML")
summary(fev1.ar1)## Generalized least squares fit by maximum likelihood
## Model: l.fev1 ~ base.age + c.age
## Data: fev1
## AIC BIC logLik
## -3549.482 -3521.492 1779.741
##
## Correlation Structure: AR(1)
## Formula: ~1 | id
## Parameter estimate(s):
## Phi
## 0.84388
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.3975811 0.05143645 -7.72956 0
## base.age 0.0962926 0.00610036 15.78475 0
## c.age 0.0867535 0.00126949 68.33736 0
##
## Correlation:
## (Intr) base.g
## base.age -0.983
## c.age -0.213 0.122
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -6.9871403 -0.3398216 0.2955394 0.9289907 3.7392820
##
## Residual standard error: 0.1682115
## Degrees of freedom: 1994 total; 1991 residualanova(fev1.ar1)## Denom. DF: 1991
## numDF F-value p-value
## (Intercept) 1 8134.326 <.0001
## base.age 1 55.843 <.0001
## c.age 1 4669.995 <.0001fev1.exch=gls(l.fev1~base.age+c.age,data=fev1,correlation = corCompSymm (form = ~ 1 | id),method="ML")
summary(fev1.exch)## Generalized least squares fit by maximum likelihood
## Model: l.fev1 ~ base.age + c.age
## Data: fev1
## AIC BIC logLik
## -2733.243 -2705.254 1371.622
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | id
## Parameter estimate(s):
## Rho
## 0.6306421
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.3202265 0.05167466 -6.19697 0
## base.age 0.0915820 0.00614303 14.90828 0
## c.age 0.0857247 0.00074700 114.75920 0
##
## Correlation:
## (Intr) base.g
## base.age -0.986
## c.age -0.124 0.071
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -7.18644465 -0.53769923 0.08346109 0.72605748 3.51756260
##
## Residual standard error: 0.1683217
## Degrees of freedom: 1994 total; 1991 residualanova(fev1.exch)## Denom. DF: 1991
## numDF F-value p-value
## (Intercept) 1 8274.091 <.0001
## base.age 1 45.867 <.0001
## c.age 1 13169.675 <.0001anova(fev1.exch,fev1.ind,fev1.ar1)## Model df AIC BIC logLik Test L.Ratio p-value
## fev1.exch 1 5 -2733.243 -2705.254 1371.6215
## fev1.ind 2 4 -1594.150 -1571.759 801.0751 1 vs 2 1141.093 <.0001
## fev1.ar1 3 5 -3549.482 -3521.492 1779.7410 2 vs 3 1957.332 <.0001The liklihood ratio statistic indiciates that the exchangeable correlation structure fits the model better than the other two correlation structures, however, the AIC and BIC were the lowest for the model with an independent correlation structure. The \(\hat{B}\)’s obtained are almost identical but the model with the exchangeable correlation structure had the smallest standard errors. The interpreation of \(\hat{B_{c}}\) is the associated effect of a one-year age difference on log(fev1) between subjects while \(\hat{B_{l}}\) is the associated effect of a one-year age difference on log(fev1) within a subject. So consistent among the models, a one-year increase in age is associated with a .092 (\(\hat{B_{c}}\)) increase in log(fev1) units across subjects. And a one-year increase in age within a subject is associated with a .086 (\(\hat{B_{l}}\)) increase in log(fev1) units.
c. Only age as a covariate
#GLS with just age as a covariate
fev1.ind.age=gls(l.fev1~age,data=fev1,method="ML")
summary(fev1.ind.age) #same as OLS as expected## Generalized least squares fit by maximum likelihood
## Model: l.fev1 ~ age
## Data: fev1
## AIC BIC logLik
## -1593.18 -1576.386 799.59
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.27415161 0.014204180 -19.30077 0
## age 0.08669269 0.001092802 79.33066 0
##
## Correlation:
## (Intr)
## age -0.967
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -7.47389430 -0.62201592 0.01233673 0.69096950 3.59045752
##
## Residual standard error: 0.162036
## Degrees of freedom: 1994 total; 1992 residualanova(fev1.ind.age)## Denom. DF: 1992
## numDF F-value p-value
## (Intercept) 1 50424.83 <.0001
## age 1 6293.35 <.0001fev1.ar1.age=gls(l.fev1~age,data=fev1,correlation = corAR1 (form = ~ 1 | id),method="ML")
summary(fev1.ar1.age)## Generalized least squares fit by maximum likelihood
## Model: l.fev1 ~ age
## Data: fev1
## AIC BIC logLik
## -3549.017 -3526.625 1778.508
##
## Correlation Structure: AR(1)
## Formula: ~1 | id
## Parameter estimate(s):
## Phi
## 0.8441701
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.3215129 0.017158098 -18.73826 0
## age 0.0869263 0.001265508 68.68884 0
##
## Correlation:
## (Intr)
## age -0.884
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -6.9213475 -0.3363911 0.2737771 0.9242987 3.7199736
##
## Residual standard error: 0.1684373
## Degrees of freedom: 1994 total; 1992 residualanova(fev1.ar1.age)## Denom. DF: 1992
## numDF F-value p-value
## (Intercept) 1 8109.537 <.0001
## age 1 4718.157 <.0001fev1.exch.age=gls(l.fev1~age,data=fev1,correlation = corCompSymm (form = ~ 1 | id),method="ML")
summary(fev1.exch.age)## Generalized least squares fit by maximum likelihood
## Model: l.fev1 ~ age
## Data: fev1
## AIC BIC logLik
## -2734.332 -2711.94 1371.166
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | id
## Parameter estimate(s):
## Rho
## 0.6316071
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -0.27228701 0.01220589 -22.30784 0
## age 0.08575976 0.00074586 114.98105 0
##
## Correlation:
## (Intr)
## age -0.732
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -7.14617793 -0.53945421 0.07136997 0.72234872 3.50165275
##
## Residual standard error: 0.168535
## Degrees of freedom: 1994 total; 1992 residualanova(fev1.exch.age)## Denom. DF: 1992
## numDF F-value p-value
## (Intercept) 1 8247.362 <.0001
## age 1 13220.643 <.0001anova(fev1.exch.age,fev1.ind.age,fev1.ar1.age)## Model df AIC BIC logLik Test L.Ratio p-value
## fev1.exch.age 1 4 -2734.332 -2711.940 1371.166
## fev1.ind.age 2 3 -1593.180 -1576.386 799.590 1 vs 2 1143.152 <.0001
## fev1.ar1.age 3 4 -3549.017 -3526.625 1778.508 2 vs 3 1957.836 <.0001The estimates among the three models were again almost identical with the exchangeavle correlation structure model having the smallest standard errors for the estimates. This model does not adjust for age at baseline/is not contstrained like the previous model, but the estimate of the association between a one-year increase in age and log(fev1) is essentially identical for these three models and the three contrained models. I would argue that the cross-sectional effect of age does not need to be adjusted for. From the scatter plots above, the mean log(fev1) is fairly constant with increases in baseline age.
d. Choosing best model
#Model diagnostics
plot(fev1.ind,id~resid(.),main="Model 1 Ind")
plot(fev1.ind.age,id~resid(.),main="Model 2 Ind")
plot(fev1.exch,id~resid(.),main="Model 1 Exch")
plot(fev1.exch.age,id~resid(.),main="Model 2 Exch")
plot(fev1.ar1,id~resid(.), main="Model 1 AR1")
plot(fev1.ar1.age,id~resid(.),"Model 2 AR1")## Warning in panel.abline(v = abl[1], ...): NAs introduced by coercion## Warning in unit(x0, default.units): NAs introduced by coercion## Warning in unit(x1, default.units): NAs introduced by coercion
plot (fev1.exch, resid (.) ~ age, abline = 0,main="Model 1 Exch")
plot (fev1.ind, resid (.) ~ age, abline = 0, main="Model 1 Ind")
plot (fev1.ar1, resid (.) ~ age, abline = 0,main="Model 1 AR1")
plot (fev1.exch.age, resid (.) ~ age, abline = 0,main="Model 2 Exch")
plot (fev1.ind.age, resid (.) ~ age, abline = 0, main= "Model 2 Ind")
plot (fev1.ar1.age, resid (.) ~ age, abline = 0,main="Model 2 AR1")
qqnorm (fev1.exch, abline = c(0,1),main="Model 1 Exch")
qqnorm (fev1.ind, abline = c(0,1), main="Model 1 Ind")
qqnorm (fev1.ar1, abline = c(0,1),main="Model 1 AR1")
qqnorm (fev1.exch.age, abline = c(0,1),main="Model 2 Exch")
qqnorm (fev1.ind.age, abline = c(0,1), main= "Model 2 Ind")
qqnorm (fev1.ar1.age, abline = c(0,1),main="Model 2 AR1")
anova(fev1.ind,fev1.ind.age,fev1.exch,fev1.exch.age,fev1.ar1,fev1.ar1.age) ## Model df AIC BIC logLik Test L.Ratio p-value
## fev1.ind 1 4 -1594.150 -1571.759 801.0751
## fev1.ind.age 2 3 -1593.180 -1576.386 799.5900 1 vs 2 2.9701 0.0848
## fev1.exch 3 5 -2733.243 -2705.254 1371.6215 2 vs 3 1144.0630 <.0001
## fev1.exch.age 4 4 -2734.332 -2711.940 1371.1660 3 vs 4 0.9110 0.3398
## fev1.ar1 5 5 -3549.482 -3521.492 1779.7410 4 vs 5 817.1499 <.0001
## fev1.ar1.age 6 4 -3549.017 -3526.625 1778.5083 5 vs 6 2.4654 0.1164Model diagnostics and an anova comparing all the models were ran. Based on the anova table, the AIC and BIC were the lowest for the model with only the age covariate that had an independence correlation structure but this has to do with the df. The AR(1) model for both sets of variables performs the worst based on AIC and BIC. I would choose the exchangeable correlation structure models because there seems to be strong within subject correlation that does not dissipate over time although the AIC/BIC is higher for than the independent structure models. Among the covariates, based on this data I would choose the unconstrained model with only age that ignores any baseline differences in log(fev1). However, given the pairwise correlations between height, age and log(fev1) it may be worthwhile to include height in the model along with age.