[R-sig-ME] level 1 variance-covariance structure

Tue Apr 12 13:25:56 CEST 2011

Dear Sebastian,

You don't need to create dummy variables your selve.

You can write m2a <- lme(attit ~ 1 +  age13 , data=data, random= ~ 0 + ind1+ ind2+ ind3+ ind4+ ind5 | id, method="REML") as

m2a <- lme(attit ~ 1 +  age13 , data=data, random= ~ 0 + factor(ind) | id, method="REML")

Or if ind is an indicator for age13:

m2a <- lme(attit ~ 1 +  age13 , data=data, random= ~ 0 + factor(age13) | id, method="REML")

Have a look at lmeControl() to increase the number of iterations. 

Best regards,

Thierry

----------------------------------------------------------------------------
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek
team Biometrie & Kwaliteitszorg
Gaverstraat 4
9500 Geraardsbergen
Belgium

Research Institute for Nature and Forest
team Biometrics & Quality Assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium

tel. + 32 54/436 185
Thierry.Onkelinx at inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

> -----Oorspronkelijk bericht-----
> Van: r-sig-mixed-models-bounces at r-project.org 
> [mailto:r-sig-mixed-models-bounces at r-project.org] Namens 
> Sebastián Daza
> Verzonden: maandag 11 april 2011 18:44
> Aan: R-SIG-Mixed-Models at r-project.org
> Onderwerp: [R-sig-ME] level 1 variance-covariance structure
> 
> Hi everyone,
> I am trying to reproduce some results models from HLM (HMLM) 
> to contrast different specifications of level 1 
> variance-covariance, but I get convergence errors. I would 
> like to know if there are any problems with my model specification...
> 
> 
> # database structure
> 
> head(data[,c(1,2,6, 9:13,17)])
>          id attit age13 ind1 ind2 ind3 ind4 ind5 ind
> 1        3  0.11    -2    1    0    0    0    0   1
> 2        3  0.20    -1    0    1    0    0    0   2
> 3        3  0.00     0    0    0    1    0    0   3
> 4        3  0.00     1    0    0    0    1    0   4
> 5        3  0.11     2    0    0    0    0    1   5
> 6        8  0.29    -2    1    0    0    0    0   1
> 
> # attit is a deviant measure and ind variables indicate 
> different waves # following some examples of snijders and 
> bosker's book, I get the unrestricted model:
> 
>  > m2a <- lme(attit ~ 1 +  age13 , data=data, random= ~ 0 + 
> ind1+ ind2+ 
> ind3+ ind4+ ind5 | id, method="REML")
> 
>  > summary(m2a)
> Linear mixed-effects model fit by REML
>   Data: data
>          AIC       BIC   logLik
>    -326.2096 -236.5348 181.1048
> 
> Random effects:
>   Formula: ~0 + ind1 + ind2 + ind3 + ind4 + ind5 | id
>   Structure: General positive-definite, Log-Cholesky parametrization
>           StdDev     Corr
> ind1     0.17219431 ind1  ind2  ind3  ind4
> ind2     0.19789253 0.493
> ind3     0.25942942 0.425 0.544
> ind4     0.28354459 0.442 0.442 0.723
> ind5     0.29097082 0.498 0.474 0.639 0.808
> Residual 0.07457025
> 
> Fixed effects: attit ~ 1 + age13
>                  Value   Std.Error  DF  t-value p-value
> (Intercept) 0.3210558 0.012832840 839 25.01829       0
> age13       0.0593529 0.004716984 839 12.58282       0
>   Correlation:
>        (Intr)
> age13 0.504
> 
> Standardized Within-Group Residuals:
>          Min          Q1         Med          Q3         Max
> -1.46371871 -0.27170442 -0.04080686  0.26239553  1.69883910
> 
> Number of Observations: 1079
> Number of Groups: 239
> 
> # variance-covariance matrix
> 
>  > extract.lme.cov2(m2a,data)$V[[6]]
>             25         26         27         28         29
> 25 0.03521160 0.01681647 0.01899029 0.02159300 0.02494013
> 26 0.01681647 0.04472218 0.02793174 0.02481343 0.02727012
> 27 0.01899029 0.02793174 0.07286434 0.05318967 0.04823107
> 28 0.02159300 0.02481343 0.05318967 0.08595826 0.06667139
> 29 0.02494013 0.02727012 0.04823107 0.06667139 0.09022474
> 
> # I get the same results than unrestricted model in HLM
> 
> # When I try to get the same unrestricted model using "corStruc" 
> commands in lme, I get a convergence problem. Am I 
> reproducing the model m2a?
> 
>  > m2b <- lme(attit ~ 1 +  age13 , data=data, random= ~ age13 
> | id, correlation = corSymm(, form =  ~ ind | id)) Error in 
> lme.formula(attit ~ 1 + age13, data = data, random = ~age13 |  :
>    nlminb problem, convergence error code = 1
>    message = iteration limit reached without convergence (9)
> 
> # When I try to get an autoregressive model,  I get again a 
> convergence problem.
> 
>  > m3a <- lme(attit ~ 1 +  age13 , data=data, random= ~ age13 
> | id, correlation = corAR1(, form =  ~ ind | id)) Error in 
> lme.formula(attit ~ 1 + age13, data = data, random = ~age13 |  :
>    nlminb problem, convergence error code = 1
>    message = iteration limit reached without convergence (9)
> 
> Does anyone know how I can solve this?
> Thank you in advance.
> 
> --
> Sebastián Daza
> sebastian.daza at gmail.com
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list 
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 