Skip to content

the equation of a mixed model fitted using nlme

2 messages · jos matejus, Thierry Onkelinx

Original
jos matejus
# 
Dear R mixed modellers

I am writing to the R mixed modelling community to request some help and
advice regarding reporting a mixed model in a publication. I have recently
received referees comments regarding a paper I submitted some time ago. The
referee has requested that I write one of the models used to analyse my
data 'statistically'. I think they mean that I should write out the
equation, and while I don?t think this is unreasonable I am having trouble
doing so. I know many of you have heard the excuse ?I am not a
statistician? but I?m afraid that this applies in my case, but I have tried
for a couple of days to figure this one out and harassed many colleagues
but without success. Therefore I am hoping I can prevail on the kindness of
the R mixed modelling community to help me with this problem.

The code used for the model was

vf1Exp<-  varExp(form=~I(days-7)|habitat*treat)

final.model <- lme(length.sq~ poly(I(days-7),2)*treat*habitat, data=mydata,
method="REML",random=~poly(I(days-7),2)|family, weights=vf1Exp)

length.sq = square root transformed length of fish

days = day following exposure (10, one day intervals starting from day 7
after exposure)

treat = two level treatment factor

habitat = two level habitat factor

family = 19 level factor



I have 10 fish per treatment combination (treat*habitat) at each time point
for each family.



The second order polynomial term for day was included to account for non
linear growth and the variance structure to account for an increase in
variance over time that was different depending on the treatment
combination. The 3-way interaction was significant.

How should I represent this model as an equation?

Thanks a million for your help.

Jos
1 day later
Thierry Onkelinx
# 
Dear Jos,

Here is an attempt. Have a look at Pinheiro & Bates (2000) for examples.

$c$ days - 7
$p1$ first order polynomial of days - 7
$p2$ second order polynomial of days - 7
$i$ family index
$j$ habitat index
$k$ treatment index
$l$ observation index

$length_{ijkl} = \beta_0 + \beta_1 p1_{ijkl} + \beta_2 p2_{ijkl} + \beta_j
+ \beta_k + \beta_{1j} p1_{ijkl}  + \beta_{2j} p2_{ijkl} + \beta_{1k}
p1_{ijkl}  + \beta_{2k} p2_{ijkl} + \beta_{1jk} p1_{ijkl}  + \beta_{2jk}
p2_{ijkl} + b_{0i} + b_{1i} p1_{ijkl} + b_{2i} p2_{ijkl} + \epsilon_{ijkl}$

$b_i = \left[{\\begin{array}{c}
  b_{0i} \\\\
  b_{1i} \\\\
  b_{1i}
\\end{array}
}\right]
\sim N(0, \Psi)$

$E[\epsilon_{ijkl}] = 0$
$Var(\epsilon_{ijkl}) = \sigma ^ 2 exp(2 \delta_{jk} c_{ijkl})$

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

2015-03-18 10:50 GMT+01:00 jos matejus <matejus106 at googlemail.com>: