Getting NaNs from nested LME

Question

Having some issues with a subset of my data (subset still has 600 values). For the experiment, I have two time points, nested within each are three treatments (TT), with 5 replicate cultures nested within them (A through E). Within each of these cultures are 20 values for individual organisms. For this subset I want to look at difference between treatments and differences between timepoints within the same treatment.

I am using R 3.1.2 and nlme package

My code is as follows:

model4a <-lme(Velocity~WeekTT, random=~1|Week/TT/Culture, method = "REML", data=body2, weights = varIdent(form=~1|WeekCulture), control=lmeControl(opt="optim"))

STR is below:

    data.frame':   600 obs. of  4 variables:
     $ Culture  : Factor w/ 5 levels "A","B","C","D",..: 1 2 3 4 5 1 2 3 4 5 
     $ Treatment: Factor w/ 3 levels "T20","T25","T25F": 1 1 1 1 1 2 2 2 2 2 
     $ Week     : num  10 10 10 10 10 10 10 10 10 10 ...
     $ Velocity : num  259 279 265 275 256 ...

Here is a screengrab of the results I get and the errors at the bottom (Working on posting some reproducible data but figured it might be a simple code error).

I've been crunching different models for the past 3 week and I think it's an easy question just my brain is frazzled and i'm overthinking it.

Please include sample input data to make a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). — MrFlick, Mar 26 '15 at 19:32

score 1 · Answer 1 · answered Mar 26 '15 at 22:30

If WeekTT is the week-by-treatment interaction, then you've included it twice in your model, once as a fixed effect and once as a random effect, because the random effects grouping structure Week/TT/Culture expands to Week + Week:TT + Week:TT:Culture, so it's not surprising that these fixed effect parameters will be unidentifiable. Perhaps, if you want to keep the fixed effect, include only the three-way interaction as the random effect.

model4a <-lme(Velocity~WeekTT, 
               random=~1|Week:TT:Culture, 
      method = "REML", data=body2, 
    weights = varIdent(form=~1|WeekCulture), 
    control=lmeControl(opt="optim"))

Depending on what the weights statement is doing (does WeekCulture represent the three-way interaction? Will a data transformation take care of the heteroscedasticity?), you might be able to just aggregate the lowest level (i.e. compute the mean of each Week:TT:Culture combination), as would be suggested by Murtaugh (2007) "Simplicity and complexity in ecological data analysis". Then you wouldn't need a mixed model at all ...

Getting NaNs from nested LME

1 Answers1