restimationgamlss

Why does gamlss give incorrect estimates of ex-gaussian distribution parameters?


From the gamlss.dist page for exGAUSS:

The ex-Gaussian distribution is often used by psychologists to model response time (RT). It is defined by adding two random variables, one from a normal distribution and the other from an exponential. The parameters mu and sigma are the mean and standard deviation from the normal distribution variable while the parameter nu is the mean of the exponential variable.

Here is how we're supposed to estimate the parameters:

library(gamlss)
y <- rexGAUS(100, mu = 300, nu = 100, sigma = 35)
m1 <- gamlss(y ~ 1, family = exGAUS)
m1

Unfortunately the estimates are way off:

Family:  c("exGAUS", "ex-Gaussian") 
Fitting method: RS() 

Call:  gamlss(formula = y ~ 1, family = exGAUS) 

Mu Coefficients:
(Intercept)  
      302.9  
Sigma Coefficients:
(Intercept)  
      3.496  
Nu Coefficients:
(Intercept)  
       4.63  

A package that has disappeared from CRAN, retimes, can still be installed from

https://cran.r-project.org/src/contrib/Archive/retimes/retimes_0.1-2.tar.gz
It has a function mexgauss:

library(retimes)
mexgauss(y)

gives:

       mu     sigma       tau 
319.42880  55.51562  85.94403 

which is closer.


Solution

  • The estimates for sigma and nu do appear to be way off in the model output. This is because exGaus() uses a log link for both those parameters by default. From the documentation, showing the defaults:

    exGAUS(mu.link = "identity", sigma.link = "log", nu.link = "log")

    The output shows results on the model scale so the estimates for those two parameters are given on the log scale.

    If we exponentiate the results for sigma and nu then we get estimates that look more reasonable and much closer to what we'd expect.

    # estimated sigma should be close to 35
    exp(3.496)
    #> [1] 32.98325
    
    # estimated nu shoud be close to 100
    exp(4.63)
    #> [1] 102.5141
    

    Created on 2021-10-11 by the reprex package (v2.0.0)