[R-sig-ME] Complex model yields similar results to simpler model, but also warnings: can I ignore them?
Thierry Onkelinx
thierry.onkelinx at inbo.be
Tue Apr 5 10:20:27 CEST 2016
Hi Jackie,
IMHO it is safer to drop the models with the warnings and keep a simpler
model. Although you have some visual evidence for the more complex models,
the data doesn't provide the evidence. This is probably because you have
too few observations for the complex model.
Best regards,
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
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
2016-03-31 18:40 GMT+02:00 Jackie Wood <jackiewood7 op gmail.com>:
> Hi Thierry,
>
> Thanks for the quick response! I've tried the following variance weight
> structures:
>
> vf1 <- varIdent (form = ~ 1 | cross*age_bin)
>
> vf2 <- varExp (form = ~ age_bin | cross)
>
> vf3 <- varComb (varIdent(form =~ 1 | cross), varExp (form = ~ age_bin))
>
> vf4 <- varExp (form = ~ age_bin)
>
> vf2 and vf3 yielded the same warnings as for the original model (using
> vf1). vf4 ran fine but I do think the data indicate that there are some
> differences in heteroscedasticity among the crosses as well. The AIC value
> was quite a bit lower for the model using vf1 but perhaps that is not
> entirely trustworthy given the warnings that were generated.
>
> I guess my original question stands about whether this warning means that
> the model is bad news and I should go with something like varExp (form = ~
> age_bin) regardless of which specification seems preferred via comparison
> of AIC values. Or if, because all other relevant output stays fairly
> consistent regardless of the variance weighting specification, I can
> proceed with the more complex version.
>
> Jackie
>
>
>
>
>
> On Thu, Mar 31, 2016 at 10:53 AM, Thierry Onkelinx <
> thierry.onkelinx op inbo.be> wrote:
>
>> Dear Jackie,
>>
>> I presume that the heteroscedasticy along age_bin is somewhat smooth. In
>> such case you use a less parametrised model the variance. Like e.g.
>> varExp(form = ~age_bin|cross). ?varClasses gives an overview of available
>> classes. Note that you can combine classes with varComb().
>>
>> Best regards,
>>
>> ir. Thierry Onkelinx
>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
>> and Forest
>> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
>> Kliniekstraat 25
>> 1070 Anderlecht
>> Belgium
>>
>> 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
>>
>> 2016-03-31 16:14 GMT+02:00 Jackie Wood <jackiewood7 op gmail.com>:
>>
>>> Hello R-users,
>>>
>>> I'm attempting to model differences in fork length over time for 4
>>> different cross types of a species of freshwater fish using the most
>>> recent
>>> version of nlme . Examining plots of the data, fork length increases
>>> non-linearly over time so I've included a second order polynomial for age
>>> such that the fixed effects portion of the model has the following
>>> specification:
>>>
>>> model <- lme(FL ~ cross * age_bin + cross*I(age_bin^2)
>>>
>>> Plots of the random effects suggest evidence for random slopes with
>>> respect
>>> to family for age and age^2, and further these are correlated with the
>>> intercept.
>>>
>>> So I specified the random effects part as:
>>>
>>> random = ~age_bin + I(age_bin^2)|fam
>>>
>>> Likelihood ratio tests do favor this random effects structure over
>>> simpler
>>> structures.
>>>
>>> Plotting the residuals, variance in length definitely increases with
>>> increasing age and also appears to vary per cross type so I added the
>>> following variance weighting structure to the model:
>>>
>>> weights = varIdent(form = ~ 1|cross*age_bin))
>>>
>>> I've performed typical likelihood ratio tests which consistently favor
>>> the
>>> model described above over other simpler model specifications (in terms
>>> of
>>> random effects specifications, and variance weighting), but with the
>>> above
>>> model I also get a few of these types of warnings:
>>>
>>> 1: In logLik.reStruct(object, conLin) :
>>> Singular precision matrix in level -1, block 15
>>>
>>> Searching online help forums the advice I see is that the model is likely
>>> overparameterized, and indeed if I remove either the variance weighting
>>> completely, or simplify the random effects to 1|fam (any random slope
>>> type
>>> random effects specification gives the same warning), everything works
>>> just
>>> fine. I also checked the raw data which seems sound to me.
>>>
>>> I do feel as though the more complex random effects structure is
>>> warranted
>>> from plotting the data and there is definitely heteroscedasticity to
>>> account for. When I run the above model without the variance weights, the
>>> resulting fixed effects coefficients and estimated random effects and
>>> correlations have values that are pretty close to the model with variance
>>> weights (the residual variance is of course different). So my question is
>>> how important are the warnings? If the output seems reasonable and
>>> corresponds pretty closely with the output of a simpler model that runs
>>> just fine, is it justifiable to ignore the warnings or am I asking for
>>> trouble?
>>>
>>> I mean, I could get rid of the variance weighting structure and simply
>>> transform fork length, it does help the heteroscedasticity issue, but I
>>> do
>>> find the variance interesting and transforming it away wouldn't be my
>>> first
>>> choice.
>>>
>>> I'd really appreciate your input!
>>>
>>> Jackie
>>>
>>> --
>>> Jacquelyn L.A. Wood, PhD.
>>> 224 Montrose Avenue
>>> Toronto, ON
>>> M6G 3G7
>>> Phone: (514) 293-7255
>>>
>>> [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-mixed-models op r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>
>>
>>
>
>
> --
> Jacquelyn L.A. Wood, PhD.
> 224 Montrose Avenue
> Toronto, ON
> M6G 3G7
> Phone: (514) 293-7255
>
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list