# [R-sig-ME] the equation of a mixed model fitted using nlme

Thierry Onkelinx thierry.onkelinx at inbo.be
Thu Mar 19 14:10:27 CET 2015

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 op googlemail.com>:

> 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
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

[[alternative HTML version deleted]]