# R gamlss - fitting and simulating lognormal response

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

• Could you tell us on what basis you expect 7000 and 36700? If you're just exponentiating the estimates for the logarithms, that's your error. See en.wikipedia.org/wiki/Log-normal_distribution for correct formulas for mean and variance in terms of the lognormal parameters.
– whuber
May 6, 2020 at 21:15
• In the step where I create df, I defined mean = 7000 and s = 36700 for x = 2. My expectation is that the using the parameters from the gamlss model and simulating random deviations from that will result in roughly the same mean and standard deviation (given some error of course). But my simmed values do not match up. I used this reference May 6, 2020 at 21:23
• How do you know your random number generation is correct? The formula you use is unusual (and I can't make any sense of it, since sig is undefined).
– whuber
May 6, 2020 at 21:30
• Sorry sig should be s. Fixed. I definitely cannot guarantee the generation is correct. Testing with sim <- rlnorm(1e6, log(7000^2 / sqrt(36700^2 + 7000^2)), sqrt(log(1 + (36700^2 / 7000^2)))) returns a mean near 7000 and sqrt(var()) near 36700 though May 6, 2020 at 21:34

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

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
$$$$
`