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  7.193461 > df_gam$sigma.fv 1.93269
Now I want to simulate y using those parameters.
sim_yhat <- (rlnorm(1e6, df_gam$mu.fv, df_gam$sigma.fv))
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)  8613.995 > sqrt(var(sim_yhat))  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.
gamlss_5.1-6 and R 4.0.0