R gamlss - fitting and simulating lognormal response

Question

I believe I am making a mistake in parametrization in this case. My goal is fit a lognormal model to data using gamlss in R, then simulate from that fitted model. Here's an example where the mean and standard deviation vary by the variable x.

library(gamlss)

#create dataset  
df <- data.frame(x = rep(0:9,1000), 
        mean = rep(c(6900,6900,7000,7600,7200,7900,7900,8100,8500,8800), 1000), 
        s = rep(c(43400,40200,36700,94200,31100,50600,45600,43600,53300,38400), 1000)) %>% 
      rowwise() %>% 
      mutate(y = rlnorm(1, log(mean^2 / sqrt(s^2 + mean^2)), sqrt(log(1 + (s^2 / mean^2))))) %>% 
      ungroup()
# fit model
df_gam <- gamlss(y ~ cs(x, df = 4), 
                 sigma.formula =~ x,
                 data = df,
                 family = LOGNO(),
                 #method = CG(),
                 trace = TRUE)

I'll use x = 2 as an example for my confusion. We get

> df_gam$mu.fv[3]
[1] 7.193461
> df_gam$sigma.fv[3]
1.93269 

Now I want to simulate y using those parameters.

sim_yhat <- (rlnorm(1e6, df_gam$mu.fv[3], df_gam$sigma.fv[3]))

But if we look at the mean and std on the response scale again, it doesn't make sense. We expect approximately 7000 and 36700, respectively, but we get

> mean(sim_yhat)
[1] 8613.995
> sqrt(var(sim_yhat))
[1] 54731.78

I believe I am missing a transformation of some sort when using the values back in rlnorm but I cannot figure out what.

Using gamlss_5.1-6 and R 4.0.0

Answer 1

I've hit this problem before as well. The documentation in GAMLSS is not really that clear about the transformations necessary to recover parameters. I've pasted an old example I made to confirm parameters are appropriately recovered. Sadly, I did not document how I landed on the transformations.

While the code below does not explicitly re-simulate the data, it should be simple to do so once the parameters are recovered correctly.

#**************************************8
## Log Normal ----
#*************************************
norm.mean=.3
norm.var=.5
lnorm.mean.true=exp(norm.mean+norm.var/2)
lnorm.var.true=exp(norm.var)*(exp(norm.var)-1)*exp(2*norm.mean)

t1=rLNO(10000,mu=norm.mean,sigma =sqrt(norm.var))

rbind(lnorm.mean.sim = mean(t1),lnorm.mean.true) # check if means match
rbind(lnorm.var.sim = var(t1), lnorm.var.true) # check if variances match

mod1=gamlss(formula = t1~1, family = LNO('identity','identity'),mu.link='identity',sigma.link='identity')
norm.mean.mod=mod1$mu.coefficients; norm.var.mod=mod1$sigma.coefficients^2

rbind(
  fit = c("mean" = norm.mean.mod,"var"= norm.var.mod),
  true = c("mean" = norm.mean[1], "var" = norm.var)
)

lnorm.mean.fit=exp(norm.mean.mod+norm.var.mod/2)
lnorm.var.fit=exp(2*norm.mean.mod+norm.var.mod)*(exp(norm.var.mod)-1)

rbind(lnorm.mean.sim = mean(t1),lnorm.mean.true, lnorm.mean.fit) # check if means match
rbind(lnorm.var.sim = var(t1), lnorm.var.true, lnorm.var.fit) # check if variances match
```