R-sig-mixed-models Digest, Vol 175, Issue 23
Not a complete answer, but some thoughts that other people might comment on. If you're primarily interested in treatment effects, because your treatments are applied at the level of fields, it might make your life a lot easier to compute average values by field & day (see Murtaugh 2007). Your data are (1) Gaussian (or being treated that way) (2) balanced. I can't see from your post whether the fixed effect model is additive (Treatment + poly(DAT, 2)) or interactive (Treatment*poly(DAT,2)). You don't really need emmeans/emtrends to test the overall effect of treatment: you could use anova() to compare the models NumberInsects ~ Treatment * poly(DAT, 2) + ... with NumberInsects ~ poly(DAT,2) + ... It's not so important if your primary interest is in Treatment, but you have at least the potential to consider variation in the temporal patterns across plots and subplots (random effect (poly(DAT, 2) | Plot/subplot) ), if you have sufficient data ... Murtaugh, Paul A. ?Simplicity and Complexity in Ecological Data Analysis.? Ecology 88, no. 1 (2007): 56?62.
On 7/27/21 1:54 PM, Inka Willms wrote:
Hello, I have a question which will probably sound very basic to most of you. However, I am insecure about my results and rather ask before presenting and publishing something I am not 100% certain about. Therefore, I would really appreciate your help! I am working on an ecotoxicological field study where a fungicide was applied onto 4 fields (with 7 subplots) whereas nothing was applied onto 4 control fields. The experiment takes place over 43 days. I checked for spatial and temporal autocorrelation. My dependent variable is for this example the number of insects. I am unfortunately not allowed to provide data, but I believe that my questions can be answered nevertheless. DAT is a continuous variable (days of the experiment) Treatment is either TRUE or FALSE I use a linear model, because the number of insects is over 2000 I constructed this model: MyModel <- glmmTMB(NumberInsects ~ Treatment poly(DAT,2) + (1|Plot/Subplot), data=dat, REML=T, family=gaussian?(link="identity")) I checked for the model fit via Dharma residual diagnostics, including temporal and spatial correlations. I answered the question for the significance of difference of number of insects in respect to treatment and control on specific days via: A <- emmeans(MyModel, ~Treatment|DAT, at =list(DAT=c(0,8,14,28,35,43))) pairs(A) However, I have issues answering the following question: * Does treatment have an effect on the number of insects, considering the complete time course of the experiment? Is it correct to use the emtrends approach here? I understand that the slope of the specified time frame are compared, which would also decouple the results from the initial numer of insects, as the slope is independent from the intercept. This is the code that I am using to answer this question: B <- emtrends(MyModel, "Treatment", var="DAT") pairs(B) It would be wrong to answer this question with the p-value of Treatment of the summary(myModel) function, right? However, I am wondering whether the 2nd polynomial of DAT is regarded in the emtrends approach. I really hope that you can help me! Kind regards, InkaMarei
________________________________
Von: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> im Auftrag von r-sig-mixed-models-request at r-project.org <r-sig-mixed-models-request at r-project.org>
Gesendet: Dienstag, 27. Juli 2021 12:00
An: r-sig-mixed-models at r-project.org <r-sig-mixed-models at r-project.org>
Betreff: R-sig-mixed-models Digest, Vol 175, Issue 22
Send R-sig-mixed-models mailing list submissions to
r-sig-mixed-models at r-project.org
To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
or, via email, send a message with subject or body 'help' to
r-sig-mixed-models-request at r-project.org
You can reach the person managing the list at
r-sig-mixed-models-owner at r-project.org
When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-mixed-models digest..."
Today's Topics:
1. Re: Multilevel equation (Brian Hudson)
----------------------------------------------------------------------
Message: 1
Date: Mon, 26 Jul 2021 22:05:19 -0400
From: Brian Hudson <bhudson.gsu at gmail.com>
To: Thierry Onkelinx <thierry.onkelinx at inbo.be>
Cc: r-sig-mixed-models <r-sig-mixed-models at r-project.org>
Subject: Re: [R-sig-ME] Multilevel equation
Message-ID:
<CANodPHM1Rg4c-B7x_whDewOcy6mm_z09T83zTmxom=yRSZozSA at mail.gmail.com>
Content-Type: text/plain; charset="utf-8"
Thierry,
Thank you! I appreciate your help and explanation - that makes sense and I
can see where my other attempts were incorrect.
A couple questions-
1) The two lines that defined eta confused me - are they supposed to be
equal to each other? I edited the equation below such that pi has the link
function and eta has the linear equation - does that work?
2) I added epsilon for the individual error
3) There was no k index (individual level) in the equations, just i and j,
so i added some indexing in the predictors.
Does this equation make sense? Any issues with what I did? Thanks again for
your (and the community's) help.
https://quicklatex.com/cache3/9f/ql_4a4eb44285f65ea0dcaf93d551c44c9f_l3.png
$$b_i\sim \mathcal{N}(0, \sigma_s^2)$$
$$b_{ij}\sim \mathcal{N}(0, \sigma_{y}^2)$$
$$\eta_{ijk} = \beta_0 + \beta_1 \textrm{X}\textsubscript{m[i]} + \beta_2
\textrm{X}\textsubscript{y[i,j]} + \beta_3
\textrm{X}\textsubscript{s[i,j,k]} + b_i + b_{ij} + \epsilon_{ijk}$$
$$\pi_{ijk} =\frac{e_{ijk}^{\eta}}{1+e_{ijk}^{\eta}}$$
$$Y_{ijk} \sim Binom(1, \pi_{ijk})$$
On Mon, Jul 19, 2021 at 1:34 PM Thierry Onkelinx <thierry.onkelinx at inbo.be>
wrote:
Dear Brian,
I'd write it as follows. In the case of a Gaussian model, you only have to
write $Y_{ijk} \sim \mathcal{N}(\eta_{ijk}, \sigma^2)$ and drop the link
function. (And you could replace \eta with \mu). Basically, Y depends on a
distribution defined by some parameters. And these parameters might need
some further definition.
$i$: state index
$j$: year index
$k$: observation index
$X_m$: state_mnthyr_pred
$X_y$: state_year_pred
$X_s$: state_pred
$$Y_{ijk} \sim Binom(\pi_{ijk})$$
$$\eta_{ijk} = \frac{\pi_{ijk}}{1- \pi_{ijk}}$$
$$\eta_{ijk} = \beta_0 + \beta_1X_m + \beta_2 X_y + \beta_3 X_s + b_i +
b_{ij}$$
$$b_i\sim \mathcal{N}(0, \sigma_s^2)$$
$$b_{ij}\sim \mathcal{N}(0, \sigma_{y}^2)$$
Best regards,
ir. Thierry Onkelinx
Statisticus / Statistician
Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND
FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx at inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be<http://www.inbo.be>
///////////////////////////////////////////////////////////////////////////////////////////
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
///////////////////////////////////////////////////////////////////////////////////////////
<https://www.inbo.be>
Op ma 19 jul. 2021 om 17:44 schreef Brian Hudson <bhudson.gsu at gmail.com>:
Hello,
I am fitting a multilevel model in `lme4` and am having trouble writing
the
equation for it. I very much appreciate any help. The formula and code is
below, but I am not sure if the equation represents the error correctly -
do I need to include error terms or is that captured by the distributions?
I am also not sure if I am representing the logit function correctly with
the indexing or functional form.
The data are comprised of US-state months nested within US-state-years and
US-states. I include predictors at each level and a varying intercept for
both state-years and states.
The formula looks like this in R:
```
as.formula(outcome ~ state_mnthyr_pred + state_year_pred + state_pred +
(1 | state) + (1 | state_year))
```
Where the outcome is dichotomous. The state months (e.g. jan-2010,
feb-2010
... jan-2013) are nested with state years and within states.
The formula I am using can be seen here:
https://quicklatex.com/cache3/e9/ql_038eeb4e4e1b0af94d3ef69fe4ff7be9_l3.png
And the LaTeX code:
$$
\begin{aligned}
\mu &=\alpha_{j[i],k[i]} +
\beta_{0}(\operatorname{state\_mnthyr\_pred})\ \\
\alpha_{j} &\sim N \left(\gamma_{0}^{\alpha} +
\gamma_{1}^{\alpha}(\operatorname{\textrm{state\_year\_pred}}),
\sigma^2_{\alpha_{j}} \right)
\text{, for \textrm{State-Year} j = 1,} \dots \text{, J} \\
\alpha_{k} &\sim N \left(\gamma_{0}^{\alpha} +
\gamma_{1}^{\alpha}(\operatorname{\textrm{state\_pred}}),
\sigma^2_{\alpha_{k}} \right)
\text{, for State k = 1,} \dots \text{, K}\\
\pi_{i} &=\frac{e_{i}^{\mu}}{1+e_{i}^{\mu}}\\
y_{i j k} \sim & \operatorname{Binom}\left(1, \pi_{i}\right)\\
\end{aligned}
$$
I really appreciate any help. Thank you.
[[alternative HTML version deleted]]
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
[[alternative HTML version deleted]]
------------------------------
Subject: Digest Footer
_______________________________________________
R-sig-mixed-models mailing list
R-sig-mixed-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
------------------------------
End of R-sig-mixed-models Digest, Vol 175, Issue 22
***************************************************
[[alternative HTML version deleted]]
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
Dr. Benjamin Bolker Professor, Mathematics & Statistics and Biology, McMaster University Director, School of Computational Science and Engineering Graduate chair, Mathematics & Statistics