2

Despite many efforts, I cannot run a linear mixed model because of the autocorrelation term. In fact, I can't manage to code for both a nested design and random slopes within it.

For example, let's imagine some monthly captures of wild rabbits in kilograms in 5 sites during 21 years:

site<- rep(rep(c("Golden Cave","Ringo's place","Damned Dam","Knockampton","Easy Fuzzy"),each=12),21) 
year <- rep(2000:2020, each=12*5)
month <- rep(seq(1,12),21*5)
rabbit_captures <- rnorm(12*21*5, 50, 10)
dataset <- as.data.frame(cbind(site,year,month,rabbit_captures))
dataset$rabbit_captures <- as.numeric(dataset$rabbit_captures)

Then is the modeling part that I can't succeed:

library(nlme)
library(MASS)

model_lme <- lme(fixed = log(rabbit_captures) ~ site,  
                 random = ~ site|year/month, 
                 correlation = corAR1(value = 0.9, form = ~ site|year/month), 
                 data = dataset, method = "ML",  
                 control = lmeControl(opt = 'optim'))

I desperately get an error that I suppose is for the "site" variable within the autocorrelation structure term:

Error in as.character.factor(X, ...) : malformed factor

I can run the model without autocorrelation:

model_lme_wo_autocorrelation <- lme(fixed = log(rabbit_captures) ~ site,  
                                    random = ~ site|year/month, 
                                    data = dataset, method = "ML", 
                                    control = lmeControl(opt = 'optim'))

And I am really interested in including autocorrelation to compare both models.

J.marine
  • 69
  • 1
  • 8
  • 1
    this wont help you much but `?corAR1` writes *"A covariate for this correlation structure must be integer valued."* – user20650 Jan 19 '20 at 12:13
  • Thank you! This time, I didn't see this line... So, how can I add an autocorrelation structure considering this particular design? Is that impossible due to the fact that we need always a numerical value within it to fit the random slopes? – J.marine Jan 19 '20 at 12:51
  • 1
    I'm not sure, sorry but perhaps is reasonable `corAR1(0.9, form = ~ 1 | year/month)` -- hopefully a stats modelling turns up! – user20650 Jan 19 '20 at 13:19
  • Actually a good idea, thanks ! With the example it went on: "Error in logLik.reStruct(object, conLin) : NA/NaN/Inf in foreign function call (arg 3)." Dataset is not consistent to run that model. With my real data, I got: "Error in lme.formula(fixed = log(CPUE + 1) ~ caleta, random = ~caleta | : problème optim, code d'erreur de convergence = 1 message = In addition: Warning message: In logLik.reStruct(object, conLin) : Singular precision matrix in level -2, block 173." I'm now pretty sure that the model would not converge. It is all right then, model without time autocorrelation is just fine. – J.marine Jan 19 '20 at 21:21
  • I consider it solved. Anyway, from a statistical point of view, I'd be really glad to know if the following model makes sense, knowing that random and autocorrelation terms differ? model_lme <- lme(fixed = log(rabbit_captures) ~ site, random = ~ site|year/month, correlation = corAR1(value = 0.9, form = ~ 1|year/month), data = dataset, method = "ML", control = lmeControl(opt = 'optim')). – J.marine Jan 19 '20 at 21:26
  • J.marine; I think it is worth asking a different question at https://stats.stackexchange.com/ ; phrased more like your comment ^^ , as in does my model make statistical sense. – user20650 Jan 19 '20 at 21:29
  • 1
    @user20650 I will ;) – J.marine Jan 19 '20 at 21:31
  • Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/206242/discussion-between-j-marine-and-user20650). – J.marine Jan 19 '20 at 21:35

0 Answers0