title: Covariance models with lme and gls
body: <p>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?</p>
  
  <pre><code># Unstructured covariance matrix
  UN &lt;- 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 &lt;- gls(y ~ ses + time, data, corr=NULL, weights=NULL, method="REML", control=lmeControl(msMaxIter = 500, msVerbose = TRUE))
  
  # Fit Random Intercept Model (RI)
  RI &lt;- 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 &lt;- 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 &lt;- 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 &lt;- 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 &lt;- 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 &lt;- 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 &lt;- 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 &lt;- 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&lt;- 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 &lt;- 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 &lt;- 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 &lt;- 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 &lt;- 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))
  </code></pre>
  
  <p>thanks</p>