Skip to content

level 1 variance-covariance structure

9 messages · Sebastián Daza, ONKELINX, Thierry, Andrew Robinson

Original
Sebastián Daza
# 
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.
ONKELINX, Thierry
# 
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
Sebastián Daza
# 
Thank you for your reply Thierry...
Increasing the number of iterations doesn't work:

m3a <- lme(attit ~ 1 +  age13 , data=data, random= ~ age13  | id, 
correlation = corAR1(, form =  ~ ind | id), control=list(maxIter=1000, 
msMaxIter=1000, niterEM=1000))

Error in lme.formula(attit ~ 1 + age13, data = data, random = ~age13 |  :
   nlminb problem, convergence error code = 1
   message = function evaluation limit reached without convergence (9)

I have attached my database. I don't know if it is a problem of my model 
or a limitation of lme function.

The best!
Sebastian.
On 4/12/2011 6:25 AM, ONKELINX, Thierry wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      ONKELINX, Thierry
    

    

      
      #
    

  


  

      
Dear Sebastian,

The model below works fine on my computer. 

m3a <- lme(attit ~ 1 +  age13 , data=dataset, random= ~ 0+factor(age13)| id, correlation = corAR1(form =  ~ age13 | id))


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

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Andrew Robinson
    

    

      
      #
    

  


  

      
Hi Sebastian,

two things to think about here ... first, the error reported below
refers to the function evaluation limit.  You need to increase that as
well as the number if iterations.  The relevant argument for control
is msMaxEval.  Second, I've had models converge after 5000
iterations.  

Best wishes

Andrew

        
On Tue, Apr 12, 2011 at 08:42:34AM -0500, Sebasti?n Daza wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Sebastián Daza
    

    

      
      #
    

  


  

      
Thierry,
I can run this model... but what does it mean?
The correlation structure that I get is:

Correlation Structure: ARMA(1,0)
  Formula: ~age13 | id
  Parameter estimate(s):
Phi1
    0

What does zero mean? I would expect get some positive number there...

  m3a <- lme(attit ~ 1 +  age13 , data, random= ~ 0 + factor(age13)| id, 
correlation = corAR1(form =  ~ age13 | id))

            

      
      
      
    

        
Linear mixed-effects model fit by REML
  Data: data
         AIC       BIC   logLik
   -324.2096 -229.5528 181.1048

Random effects:
  Formula: ~0 + factor(age13) | id
  Structure: General positive-definite, Log-Cholesky parametrization
                 StdDev     Corr
factor(age13)-2 0.17219431 f(13)-2 f(13)-1 f(13)0 f(13)1
factor(age13)-1 0.19789254 0.493
factor(age13)0  0.25942941 0.425   0.544
factor(age13)1  0.28354459 0.442   0.442   0.723
factor(age13)2  0.29097081 0.498   0.474   0.639  0.808
Residual        0.07457025

Correlation Structure: ARMA(1,0)
  Formula: ~age13 | id
  Parameter estimate(s):
Phi1
    0
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.46371874 -0.27170442 -0.04080688  0.26239551  1.69883907

Number of Observations: 1079
Number of Groups: 239

        
On 4/12/2011 10:21 AM, ONKELINX, Thierry wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Sebastián Daza
    

    

      
      #
    

  


  

      
Thank you Andrew,
But it doesn't work, I get the same error:

  m3a <- lme(attit ~ 1 +  age13 , data=data, random= ~ age13  | id, 
correlation = corAR1(, form =  ~ ind | id), 
control=list(msMaxEval=10000, maxIter=10000, msMaxIter=10000, 
niterEM=10000))
Error in lme.formula(attit ~ 1 + age13, data = data, random = ~age13 |  :
   nlminb problem, convergence error code = 1
   message = singular convergence (7)

        
On 4/12/2011 2:05 PM, Andrew Robinson wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      ONKELINX, Thierry
    

    

      
      #
    

  


  

      
There is no auto-correlation left AFTER the fixed and random effects are taken into account. So you probably will have to choose between the models below.

m3a <- lme(attit ~ age13 , data, random= ~ 0 + factor(age13)| id)
m3b <- lme(attit ~ age13 , data, random= ~ 1| id, correlation = corAR1(form =  ~ age13 | id))

----------------------------------------------------------------------------
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

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Sebastián Daza
    

    

      
      #
    

  


  

      
Finally, I had to run the models using HLM, where I could get results 
(without iteration errors) for a model like this:

attit ~ age13 , data, random= ~ age13 | id, correlation = corAR1(form = 
  ~ wave | id))

I don't know why I could run this kind of models in HLM and I couldn't 
do it using R (lme). It would be good to know if it is a limitation of 
the R package or an over-parametrization of the model... or it is 
related to a specific characteristic of the dataset. I don't have any clue.

Thank you!

        
On 4/13/2011 3:29 AM, ONKELINX, Thierry wrote: