cats.uni <- read.csv("http://www.stat.ufl.edu/~winner/data/cat_zyl_uni.csv") attach(cats.uni); names(cats.uni) id <- factor(id) trt <- factor(trt) timepnt1 <- factor(timepnt,ordered=T) timepnt <- factor(timepnt, levels=1:5, labels=c("T1","T2","T3","T4","T5")) library(nlme) ## Simple Covariance Structure (Measurements within cats uncorrelated) fit.sim <- gls(y ~ trt*timepnt) summary(fit.sim) anova(fit.sim) AIC(fit.sim) # getVarCov(fit.sim, individual="1",type="conditional") ## Compound Symmetry estimates the 2 Variances but seperate for ZGZ' and R fit.cs <- lme(y ~ trt*timepnt, random= ~ 1 | id, corr = corCompSymm(, form= ~ 1 | id)) summary(fit.cs) anova(fit.cs) AIC(fit.cs) getVarCov(fit.cs, individual="1",type="conditional") ## Compound Symmetry that prints marginal V-Cov matrix = ZGZ'+R w/in cats fit.cs1 <- lme(y ~ trt*timepnt, random= list(id = pdCompSymm(~ trt - 1))) summary(fit.cs1) anova(fit.cs1) AIC(fit.cs1) getVarCov(fit.cs1, individual="1",type="marginal") ## Unstructured Covariance Structure fit.un <- lme(y ~ trt*timepnt, random= ~1 | id, corr = corSymm(form= ~ 1 | id), weights = varIdent(form = ~ 1 | timepnt), control=lmeControl(maxIter=100, opt="optim")) summary(fit.un) anova(fit.un) AIC(fit.un) getVarCov(fit.un, individual="1",type="conditional") ## AR(1) Covariance Structure fit.ar1 <- lme(y ~ trt*timepnt, random= ~1 | id, corr = corAR1(, form= ~ 1 | id)) summary(fit.ar1) anova(fit.ar1) AIC(fit.ar1) getVarCov(fit.ar1, individual="1",type="conditional") ## Heterogeneous AR(1) Covariance Structure fit.arh1 <- lme(y ~ trt*timepnt, random= ~1 | id, corr = corAR1(, form= ~ 1 | id), weight = varIdent(form = ~ 1 | timepnt), control=lmeControl(maxIter=100, opt="optim")) summary(fit.arh1) anova(fit.arh1) AIC(fit.arh1) getVarCov(fit.arh1, individual="1",type="conditional") ### Using GLS Covariance Structure (No Random Effects) ## Compound Symmetry fit.cs.g <- gls(y ~ trt*timepnt, # random= ~ 1 | id, corr = corCompSymm(, form= ~ 1 | id)) summary(fit.cs.g) anova(fit.cs.g) AIC(fit.cs.g) getVarCov(fit.cs.g, individual="1",type="conditional") ## Unstructured Covariance Structure fit.un.g <- gls(y ~ trt*timepnt, # random= ~1 | id, corr = corSymm(form= ~ 1 | id), weights = varIdent(form = ~ 1 | timepnt), control=lmeControl(maxIter=100, opt="optim")) summary(fit.un.g) anova(fit.un.g) AIC(fit.un.g) getVarCov(fit.un.g, individual="1",type="conditional") ## AR(1) Covariance Structure fit.ar1.g <- gls(y ~ trt*timepnt, # random= ~1 | id, corr = corAR1(, form= ~ 1 | id)) summary(fit.ar1.g) anova(fit.ar1.g) AIC(fit.ar1.g) getVarCov(fit.ar1.g, individual="1",type="conditional") ## Heterogeneous AR(1) Covariance Structure fit.arh1.g <- gls(y ~ trt*timepnt, # random= ~1 | id, corr = corAR1(, form= ~ 1 | id), weight = varIdent(form = ~ 1 | timepnt), control=lmeControl(maxIter=100, opt="optim")) summary(fit.arh1.g) anova(fit.arh1.g) AIC(fit.arh1.g) getVarCov(fit.arh1.g, individual="1",type="conditional")