lme makes Rstudio crash and pc nonresponsive after increasing memory limit

Question

I am working with a data set of 205 observations and 2 explanatory variables : site (two levels) and strain ( 21 levels ). I am trying to fit a mixed model to the data when strain is the fixed variable. the experiment is not balanced ( i.e. there are different numbers of repetitions in each strain*site combination). And I use Rstudio 3.0.2; windows 8.

I first tried to run the following -

 lme(dist_OFT_min1~strain , 
     random=~strain|site,data=data.frame(d1),
     method='REML')

Yet I got the following error :

Error in logLik.lmeStructInt(lmeSt, lmePars) : 'Calloc' could not allocate memory (730512784 of 8 bytes)
In addition: Warning messages:
1: In logLik.lmeStructInt(lmeSt, lmePars) : Reached total allocation of 3873Mb: see help(memory.size)
2: In logLik.lmeStructInt(lmeSt, lmePars) : Reached total allocation of 3873Mb: see help(memory.size)

I tried to increase the memory limit to 8000 . But then when I rerun the previous line Rstudio (and the pc in total ) get very slow and unresponsive, I let it run for over 30 minutes for no result.

I tried to fit the model on a cut-down data set - containing only 5 of the strains :

 lme(dist_OFT_min1~strain , 
     random=~strain|site,
     data=data.frame(d.test),
     method='REML',
     control = lmeControl(msMaxIter =90,maxIter=90)) 

and another kind of trouble popped :

Error in lme.formula(dist_OFT_min1 ~ strain, random = ~strain | site, :
nlminb problem, convergence error code = 1 message = iteration limit reached without convergence (10)

Can anyone please help me understand what kind of trouble I might be facing with my data ? How can I check for any ? (I'm still a beginner). Here I share similar data I simulated on my own, Yet I faced the same trouble as with the original data ( lab equals site ):

    library(nlme)
    u.strain=letters[1:10]
    u.lab=paste('L',letters[1:5],sep='')
    test.dat=data.frame(strain=sort(rep(u.strain,5)),lab=rep(u.lab,10),test=rep(7,50))
    test.dat[order(test.dat$strain),'test']=test.dat[order(test.dat$strain),'test']+sort(rep(rnorm(n=10,mean=0,sd=3),5)) #strain effect
    test.dat[order(test.dat$lab),'test']=test.dat[order(test.dat$lab),'test']+sort(rep(rnorm(n=5,mean=0,sd=4),10)) #lab effect
    test.dat[,'test']=test.dat[,'test']+ rnorm(n=50,mean=0,sd=5) # interaction
    test.dat=rbind(test.dat,test.dat,test.dat,test.dat,test.dat)
    test.dat[,'test']=test.dat[,'test']+rnorm(n=250,0,sd=2.5)

I figured that the model I actually need is as follows :

    lme(teat~strain , random=~strain+strain:lab|lab,
             data=data.frame(test.dat),method='REML')

Still the same trouble popped.

p.s. I mainly seek to estimate the variance of the interaction strain:lab. Any tips ?

Answer 1

Define interaction explicitly:

test.dat$strainlab <- with(test.dat,interaction(strain,lab))

Fit site (lab) as fixed effect (makes much more sense if you have only two sites), strain-by-site interaction as random:

m1 <- lme(test~strain+lab , random=~1|strainlab,
    data=data.frame(test.dat),method='REML')
VarCorr(m1)
## strainlab = pdLogChol(1) 
##             Variance  StdDev
## (Intercept) 29.286446 5.411695
## Residual     5.398621 2.323493

Approximate confidence intervals:

intervals(m1,which="var-cov")
## Approximate 95% confidence intervals
## 
##  Random Effects:
##   Level: strainlab 
##                    lower     est.    upper
## sd((Intercept)) 4.258995 5.411695 6.876374
## 
##  Within-group standard error:
##    lower     est.    upper 
## 2.106594 2.323493 2.562725 
##

Alternatively, fit lab and strain-by-lab interaction as random:

m2 <- lme(test~strain , random=~1|lab/strain,
data=data.frame(test.dat),method='REML')
VarCorr(m2)
##             Variance     StdDev  
## lab =       pdLogChol(1)         
## (Intercept) 18.608568    4.313765
## strain =    pdLogChol(1)         
## (Intercept) 29.286446    5.411695
## Residual     5.398621    2.323493
## 
intervals(m2,which="var-cov")
## Approximate 95% confidence intervals
## 
##  Random Effects:
##   Level: lab 
##                    lower     est.    upper
## sd((Intercept)) 1.925088 4.313765 9.666346
##   Level: strain 
##                    lower     est.    upper
## sd((Intercept)) 4.259026 5.411695 6.876324
## 
##  Within-group standard error:
##    lower     est.    upper 
## 2.106600 2.323493 2.562717