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