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