2
$\begingroup$

I know that there are many related questions and answers, but each setup, including my case, is so specific in LMM that it is hard to draw reliable analogies.

The data (available here):

  • subject_id: Participant ID.
  • rt_start: Response time (RT) to stimulus (between 150 and 1000 ms).
  • stim_type: Either of two types of the displayed stimulus, probe or control.
  • block_number: Block (1, 2, 3, or 4). (Different blocks contain different specific stimulus words, but their types are always either probe or control.)
  • trial_number: Trial, numbered 1-648, with ca. 46% "missing" trials randomly throughout, but relatively evenly distributed (too slow and incorrect responses, and other, irrelevant types excluded from the data, originally mixed together with probe and control stimuli in a balanced sequence).

It is important that the stim_type effect is large and expected: probe RT is on average much larger than control RT. What I want to test is just how much this difference is influenced by learning (practice/habituation effect) throughout the task, expecting that the difference is decreasing with time. Crucially, I want to see both (a) whether there is a within-block effect (increasing trial number leads to decreasing probe-control difference; stim_type x trial_number interaction), and (b) whether there is a between-block effect (increasing block number leads to decreasing probe-control difference; stim_type x block_number interaction)

I'm pretty sure that I need LMM to do this, but I'm not sure how exactly.

Some general relevant points:

  • I expect an overall learning effect (generally faster responses), hence both block_number and trial_number main fixed effects are included, but otherwise they are not variables of interest in themselves.
  • Since each trial happens in each block (in principle), perhaps one could argue that trials in the first blocks (1-162) are no different from e.g. those in the fourth blocks (1-162 in that block, but counting with the entire test 487-648). However, in my present modelling this is not taken into account.
  • I assume that the degree of change per trial and per block are both identical, hence I treat them as numeric (continuous) variables.
  • The general RT as well as the item_type factor (probe vs. control) differs by subject, hence I include (item_type|subjects) that should, if I understand correctly, account for both general RT baseline (intercept) and for probe-control difference ("slope" of the dummy variable item_type) per individual.
  • I'm using glmer with Gamma(link = "identity"), as this is recommended for RTs, plus this circumvents the normality assumption (which would be violated). I had no success with robust alternatives (e.g. robustlmm and similar), mostly because they run into errors due to the large sample.
  • No matter how I vary the approach (e.g. any of the parameters mentioned above, or even with block-wise aggregation using simple ANOVA), the overall length effect (e.g., stim_type x trial_number interaction) seems robust (e.g. at least p < .015 for everything I tried so far), so there seems to be little doubt about it. I even tried robustlmm:rlmer with decimated data (to avoid the large-data error), and even so the effects were significant and with similar estimates. What's more, I got the same results (again very low p values) with another very similar independent dataset (omitted here for brevity), where residuals look somewhat better.

library('lme4')

# fitting full model
mlm_full = glmer(
rt_start ~ stim_type + block_number + trial_number +
    stim_type:block_number + stim_type:trial_number +
    (stim_type | subject_id),
data = lgcit_dat,
family = Gamma(link = "identity")
)


> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
> Family: Gamma  ( identity )
> Formula: rt_start ~ stim_type + block_number + trial_number + stim_type:block_number +  
>     stim_type:trial_number + (stim_type | subject_id)
> Data: lgcit_dat
>     AIC       BIC    logLik  deviance  df.resid 
> 910966.3  911058.7 -455473.1  910946.3     76637 
> Random effects:
> Groups     Name           Std.Dev. Corr 
> subject_id (Intercept)    20.1051       
>             stim_typeprobe 28.2413  -0.07
> Residual                   0.1904       
> Number of obs: 76647, groups:  subject_id, 219
> Fixed Effects:
>                 (Intercept)               stim_typeprobe                 block_number  
>                 514.31118                    115.23858                     18.92389  
>             trial_number  stim_typeprobe:block_number  stim_typeprobe:trial_number  
>                 -0.19776                     -3.21046                     -0.05767  


# without stim_type:trial_number
mlm_xtrial = glmer(
rt_start ~ stim_type + block_number + trial_number +
    stim_type:block_number +
    (stim_type | subject_id),
data = lgcit_dat,
family = Gamma(link = "identity")
)

# without stim_type:block_number
mlm_xblock = glmer(
rt_start ~ stim_type + block_number + trial_number +
    stim_type:trial_number +
    (stim_type | subject_id),
data = lgcit_dat,
family = Gamma(link = "identity")
)

# test whether stim_type:trial_number is significant contributor
aov_trials = anova(mlm_full, mlm_xtrial, )

>         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
> mlm_xtrial    9 910973 911056 -455477   910955                        
> mlm_full     10 910966 911059 -455473   910946 8.3873  1   0.003779 **



# test whether stim_type:block_number is significant contributor
aov_trials = anova(mlm_full, mlm_xblock, )

>         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
> mlm_xblock    9 910965 911048 -455474   910947                     
> mlm_full     10 910966 911059 -455473   910946 0.9347  1     0.3336

Questions:

  1. Is this approach generally correct? In particular: is this how I should include and test the significance of the two critical interactions?

  2. Looking at the full model, would it be correct to say that the probe-control differences decrease significantly by 0.06 ms per each trial, and nominally (but without statistical significance) by an additional 3 ms per each block? (Note that these coefficients too hardly change depending on settings.)

$\endgroup$
2
  • 2
    $\begingroup$ If you're fitting a gamma GLM, why are you checking for residual normality? $\endgroup$ May 22, 2021 at 6:22
  • 1
    $\begingroup$ Because I'm bad at statistics ;) But seriously, good point, noted. Thanks! $\endgroup$
    – gaspar
    May 22, 2021 at 10:12

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.