# [R-sig-ME] R-sig-mixed-models Digest, Vol 175, Issue 23

Ben Bolker bbo|ker @end|ng |rom gm@||@com
Wed Jul 28 01:39:26 CEST 2021

   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 using r-project.org> im Auftrag von r-sig-mixed-models-request using r-project.org <r-sig-mixed-models-request using r-project.org>
> Gesendet: Dienstag, 27. Juli 2021 12:00
> An: r-sig-mixed-models using r-project.org <r-sig-mixed-models using 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 using 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 using r-project.org
>
> You can reach the person managing the list at
>          r-sig-mixed-models-owner using r-project.org
>
> 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 using gmail.com>
> To: Thierry Onkelinx <thierry.onkelinx using inbo.be>
> Cc: r-sig-mixed-models <r-sig-mixed-models using r-project.org>
> Subject: Re: [R-sig-ME] Multilevel equation
> Message-ID:
>          <CANodPHM1Rg4c-B7x_whDewOcy6mm_z09T83zTmxom=yRSZozSA using 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 using 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 using 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 using 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 using 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 using 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 using 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