# Analyzing reaction time data by implementing GAMM analyses with non-normal distribution parameters

I have two questions, a conceptual one and a practical one (that are closely related). And just as a note, I'm not super familiar with this level of stats (much more comfortable with simple linear regression and ANOVA), so apologies if there are oversimplifications/obvious conceptual issues at play here.

First: There are clearly issues with using means to analyze reaction time data because RT data are not normally distributed. An appropriate distribution that's frequently referenced is the ex-Gaussian distribution. Because I wanted to implement an analysis using an appropriate distribution without transforming the data I looked into using the GAMLSS package to run the analysis. But then I go back to the literature and what I see is that the ex-Gaussian parameters (mu, sigma, and tau) are being extracted from the gam model and then being used in a GLMM. So the question is: why not just report and interpret the gam model with an ex-Gaussian distribution? Is there anything incorrect about doing so in order to examine these data with a non-normal distribution?

Second: a brief rundown of the study then the questions (sample output is at the end). Participants completed a simple RT task at two timepoints at two different locations. The RT task has a condition called valence with three levels. Participants are randomized to watch 1 of 2 videos at both timepoints. Participants complete a measure called STAI once. Variables: location (w/in factor with two levels); valence (w/in factor with three levels); video (between factor with two levels); STAI (continuous variable). Random intercepts are included for individual participants.

There are a few things I don't understand about GAMLSS: 1) The output is different when I include different estimation parameters (so if I include only mu in the formula e.g. gamlss(formula = RT~location, data=data) it's different from gamlss(RT~location, sigma.fo = ~location, nu.fo = ~location, data=data) Any thoughts on why this is the case? Are fit statistics reliable between these two models? I post an example output with the mu and nu formulas only at the end. 2) This is a continuation of the first part: my understanding of interpreting the parameters is that they each refer to a distinct component (mu refers to mean; sigma to sd; tau to the exponential component of the mean and sd). So if location (which strongly influences RT in every model I try and run except when I included the tau component in the analysis) is significant for tau, but not mu - does this mean that the real explanation for the difference of location on RTs is not global (i.e. average) but somehow explained by the exponential convolution of the mean and sd? How would I explain this?

Below I provide an example output of the model I ran first with all formulas, then only the mu and nu formula (gamlss refers to nu, but it's actually equivalent to tau in the context of the exGAUS distribution, I think?).

Family:  c("exGAUS", "ex-Gaussian")

Call:
gamlss(formula = RT ~ STAI_T_tot + location + video + valence +
random(subID), sigma.formula = ~STAI_T_tot + location + video +
valence + random(subID), nu.formula = ~STAI_T_tot + location +
video + valence + random(subID), family = exGAUS(), data = na.omit(data))

Fitting method: RS()

-------------------------------------------------------------------
Mu Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)    606.19898    2.48220 244.219  < 2e-16 ***
STAI_T_tot      -0.78644    0.05277 -14.905  < 2e-16 ***
location1        1.03000    1.11498   0.924   0.3556
video1           8.08323    1.11894   7.224 5.09e-13 ***
valenceneutral   3.37995    1.35687   2.491   0.0127 *
valencepos       2.82100    1.35826   2.077   0.0378 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-------------------------------------------------------------------
Sigma Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)     4.6897578  0.0293722 159.667  < 2e-16 ***
STAI_T_tot     -0.0015505  0.0006102  -2.541  0.01106 *
location1       0.0236655  0.0131948   1.794  0.07289 .
video1          0.0384450  0.0131913   2.914  0.00356 **
valenceneutral  0.0185930  0.0161728   1.150  0.25029
valencepos      0.0255943  0.0161156   1.588  0.11225
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-------------------------------------------------------------------
Nu Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)     5.6664104  0.0150130 377.432   <2e-16 ***
STAI_T_tot     -0.0030183  0.0003188  -9.467   <2e-16 ***
location1       0.1074412  0.0070122  15.322   <2e-16 ***
video1          0.0640823  0.0070783   9.053   <2e-16 ***
valenceneutral -0.0107050  0.0085951  -1.245    0.213
valencepos     -0.0102235  0.0085242  -1.199    0.230
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-------------------------------------------------------------------
NOTE: Additive smoothing terms exist in the formulas:
i) Std. Error for smoothers are for the linear effect only.
ii) Std. Error for the linear terms may not be reliable.
-------------------------------------------------------------------
No. of observations in the fit:  80795
Degrees of Freedom for the fit:  436.9502
Residual Deg. of Freedom:  80358.05
at cycle:  20

Global Deviance:     1115911
AIC:     1116784
SBC:     1120848

Family:  c("exGAUS", "ex-Gaussian")

Call:
gamlss(formula = RT ~ STAI_T_tot + location + video + valence +
random(subID), nu.formula = ~STAI_T_tot + location + video +
valence + random(subID), family = exGAUS(), data = na.omit(data))

Fitting method: RS()

-------------------------------------------------------------------
Mu Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)    540.62127    2.58243 209.346  < 2e-16 ***
STAI_T_tot      -0.15510    0.05369  -2.889  0.00387 **
location1       -1.34292    1.17034  -1.147  0.25119
video1         -15.72136    1.17199 -13.414  < 2e-16 ***
valenceneutral   1.37373    1.43494   0.957  0.33840
valencepos       0.60667    1.43042   0.424  0.67148
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-------------------------------------------------------------------
Sigma Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.418776   0.006184   714.6   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-------------------------------------------------------------------
Nu Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)     5.8552886  0.0152149 384.838   <2e-16 ***
STAI_T_tot     -0.0045268  0.0003249 -13.934   <2e-16 ***
location1       0.1100487  0.0071404  15.412   <2e-16 ***
video1          0.0889630  0.0072172  12.327   <2e-16 ***
valenceneutral -0.0034509  0.0087375  -0.395    0.693
valencepos     -0.0029106  0.0086685  -0.336    0.737
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-------------------------------------------------------------------
NOTE: Additive smoothing terms exist in the formulas:
i) Std. Error for smoothers are for the linear effect only.
ii) Std. Error for the linear terms may not be reliable.
-------------------------------------------------------------------
No. of observations in the fit:  80795
Degrees of Freedom for the fit:  294.7012
Residual Deg. of Freedom:  80500.3
at cycle:  20

Global Deviance:     1119469
AIC:     1120059
SBC:     1122799


That works well because in the RT literature using the ex-Gaussian as the model of cognition, there are theoretical reasons to believe that experimental manipulations will affect one or more parameters by not the others (but see Matzke & Wagenmakers, 2009). In your case, one might theorize that, say, valence only affects $\mu$, but not the other parameters (I am just picking at random here). So then you'd specify an equation only for $\mu$ but not the other variables. Everything else would go in the first formula call, as you have done.