[R-sig-ME] Complex model yields similar results to simpler model, but also warnings: can I ignore them?

[Jackie Wood]

Hi Thierry,

Ok! Thanks for the advice; I'll go ahead with the simpler model.

Cheers,
Jackie


On Tue, Apr 5, 2016 at 4:20 AM, Thierry Onkelinx <thierry.onkelinx at inbo.be>
wrote:

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


-- 
Jacquelyn L.A. Wood, PhD.
224 Montrose Avenue
Toronto, ON
M6G 3G7
Phone: (514) 293-7255

	[[alternative HTML version deleted]]