Covariance models with lme and gls

Question

I am trying to find several covariance models using gls and lme. My goal is to determine which covariance model suits my data better. I am afraid, however, that I am not specifying the code correctly. Can someone take a look at my code and help me figure out if I am doing everything right 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