I'm trying to fit several covariance models using gls and lme. My aim is to
identify which covariance model fits my data better. I'm afraid, however,
that I'm not specifying the code properly. Could someone take a look on my
code and help me figuring out whether I'm pursuing everything correctly?

# Unstructured covariance matrix
UN <- gls(y ~ ses + time, data, corr=corSymm(form=~1|id),
weights=varIdent(form=~1|time), method="REML", control=lmeControl(msMaxIter
= 500, msVerbose = TRUE), na.action=na.omit)

# Independence covariance matrix
IN <- gls(y ~ ses + time, data, corr=NULL, weights=NULL, method="REML",
control=lmeControl(msMaxIter = 500, msVerbose = TRUE))

# Fit Random Intercept Model (RI)
RI <- lme(y ~ ses + time, data, na.action=na.omit, method="REML",
random=~1|id, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Fit Random Intercept and Slopes Model (RIAS)
RIAS <- lme(y ~ ses + time, data, na.action=na.omit, method="REML",
random=~time | id, control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# Fit Compound Symmetry Error Covariance Matrix
CS <- gls(y ~ ses + time, data, na.action=na.omit, method="REML",
correlation=corCompSymm(,form=~1|id), control=lmeControl(msMaxIter = 500,
msVerbose = TRUE))

# Fit Heterogeneous Compound Symmetry Error Covariance Matrix
CSH <- gls(y ~ ses + time, data, na.action=na.omit, method="REML",
correlation=corCompSymm(,form=~1|id), weights=varIdent(form=~1|time),
control=lmeControl(msMaxIter = 500, msVerbose = TRUE))

# AR(1)
AR1 <- gls(y ~ ses + time, data, na.action=na.omit, method="REML",
correlation=corAR1(,form=~1|id), control=lmeControl(msMaxIter = 200,
msVerbose = TRUE))

# AR(1) under heterocedasticity
ARH1 <- gls(y ~ ses + time, data, na.action=na.omit, method="REML",
correlation=corAR1(,form=~1|id), weights=varIdent(form=~1|time),
control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# RI plus AR(1)
RIAR1 <- lme(y ~ ses + time, data, na.action=na.omit, method="REML",
random=~1|id, correlation=corAR1(form=~1|id), control=lmeControl(msMaxIter
= 200, msVerbose = TRUE))

# RI plus AR(1) under heterocedasticity
RIARH1 <- lme(y ~ ses + time, data, na.action=na.omit, method="REML",
random=~1|id, correlation=corAR1(form=~1|id),
weights=varIdent(form=~1|time), control=lmeControl(msMaxIter = 200,
msVerbose = TRUE))

# RIAS plus AR(1)
RIASAR1<- lme(y ~ ses + time, data, na.action=na.omit, method="REML",
random=~time|id, correlation=corAR1(form=~1|id),
control=lmeControl(msMaxIter = 200, msVerbose = TRUE))

# ARMA(1,1)
ARMA11 <- gls(y ~ ses + time, data, na.action=na.omit, method="REML",
correlation=corARMA(,form=~time|id, p=1, q=1), control=lmeControl(msMaxIter
= 200, msVerbose = TRUE))

# ARMA(1,1) under heterocedasticity
ARMA11HE <- gls(y ~ ses + time, data, na.action=na.omit, method="REML",
correlation=corARMA(,form=~time|id, p=1, q=1),
weights=varIdent(form=~1|time),  control = lmeControl(msMaxIter = 200,
msVerbose = TRUE))

# Fit Toeplitz Error Covariance Matrix
TOEP <- gls(y ~ ses + time, data, na.action=na.omit, method="REML",
correlation=corARMA(,form=~1|id, p=3, q=0))

# RIAQS plus AR(1) allow for heterocedasticity
RIAQSAR1 <- lme(y ~ ses + time, data, na.action=na.omit, method="REML",
weights=varPower(form=~time), random=~time + I(time^2)|id,
correlation=corAR1(form=~time), control=lmeControl(msMaxIter = 200,
msVerbose = TRUE))

thanks

	[[alternative HTML version deleted]]


