R/est_param_2.R
est_param_2.RdThis function estimates parameters like mean, standard deviation for a gamlss2 object. This function is very useful when the mean or standard deviation do not match with the parameters \(\mu\), \(\sigma\), \(\nu\) and \(\tau\), ..., or some function of the parameters \(\mu\), \(\sigma\), \(\nu\) and \(\tau\).
est_param_2(mod, fun = "mean", m = 100, data, newdata = NULL, ...)a gamlss2 object.
the estimated parameter, by default is mean but can be changed for sd, var or any statistical function.
a number to indicate the number of observations to simulate, by default its values is 100.
the original data frame used to fit the mod.
a data frame containing new values for the explanatory variables.
additional parameters for "fun" function.
est_param_2 function returns a vector.
The function obtains the fitted values for mu, sigma, nu and tau. The functions simulates m observations and then it obtains the mean/fun defined by the user.
# Example for NOrmal regression ------------------------
# Here E(y) = mu and Var(y)= sigma
n <- 7
x <- runif(n=n)
y <- rnorm(n=n, mean=-5 + 10 * x, sd=exp(1 + 2 * x))
data1 <- data.frame(y, x)
data1
#> y x
#> 1 -0.4271114 0.05998003
#> 2 2.6699518 0.13758607
#> 3 5.1379106 0.14374278
#> 4 21.1474481 0.96052288
#> 5 13.0618165 0.56387008
#> 6 -5.6565709 0.11247740
#> 7 -15.3285801 0.80137135
library(gamlss2)
#> Warning: package ‘gamlss2’ was built under R version 4.4.3
#> Loading required package: mgcv
#> This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.
#>
#> Attaching package: ‘gamlss2’
#> The following objects are masked from ‘package:gamlss’:
#>
#> GAIC, Rsq
mod1 <- gamlss2(y ~ x, sigma.fo=~x, family=NO, data=data1)
#> GAMLSS-RS iteration 1: Global Deviance = 46.8837 eps = 0.125365
#> GAMLSS-RS iteration 2: Global Deviance = 46.829 eps = 0.001166
#> GAMLSS-RS iteration 3: Global Deviance = 46.8287 eps = 0.000006
# Obtaining the means
m1 <- fitted(mod1, model="mu", type="parameter")
m2 <- est_param_2(mod1, m=10000, fun="mean")
cbind(m1, m2) # Comparing the estimated means
#> m1 m2
#> 1 -0.4153369 -0.4288510
#> 2 0.5489152 0.5509515
#> 3 0.6254121 0.6086742
#> 4 10.7738732 11.0339829
#> 5 5.8454776 5.8443838
#> 6 0.2369410 0.2363658
#> 7 8.7964217 8.4613334
cor(m1, m2)
#> [1] 0.9993113
# Obtaining the standard deviations
s1 <- fitted(mod1, model="sigma", type="parameter")
s2 <- est_param_2(mod1, m=10000, fun="sd")
cbind(s1, s2) # Comparing the estimated standard deviations
#> s1 s2
#> 1 3.349234 3.364794
#> 2 3.950544 3.909309
#> 3 4.002635 4.027919
#> 4 22.755573 22.743399
#> 5 9.785041 9.811429
#> 6 3.745031 3.771072
#> 7 16.219019 16.090224
cor(s1, s2)
#> [1] 0.999978
# Example for GAmma regression --------------------------
# Here E(y) = mu and Sd(y) = mu * sigma
n <- 7
x <- runif(n=n)
y <- rGA(n=n, mu=exp(-1 + 5 * x), sigma=exp(-1 + 2 * x))
data2 <- data.frame(y, x)
data2
#> y x
#> 1 0.1799952 0.009674706
#> 2 0.8461014 0.764305320
#> 3 3.8110424 0.345629672
#> 4 4.8053332 0.414382810
#> 5 1.3435752 0.172801998
#> 6 1.4552517 0.827368697
#> 7 3.3387025 0.453946590
mod2 <- gamlss2(y ~ x, sigma.fo=~x, family=GA, data=data2)
#> GAMLSS-RS iteration 1: Global Deviance = 23.9863 eps = 0.056003
#> GAMLSS-RS iteration 2: Global Deviance = 20.9131 eps = 0.128125
#> GAMLSS-RS iteration 3: Global Deviance = 19.7666 eps = 0.054822
#> GAMLSS-RS iteration 4: Global Deviance = 19.7481 eps = 0.000936
#> GAMLSS-RS iteration 5: Global Deviance = 19.748 eps = 0.000002
summary(mod2)
#> Call:
#> gamlss2(formula = y ~ x, data = data2, family = GA, ... = pairlist(sigma.fo = ~x))
#> ---
#> Family: GA
#> Link function: mu = log, sigma = log
#> *--------
#> Parameter: mu
#> ---
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.9609 0.6396 3.066 0.0547 .
#> x -2.1358 0.8885 -2.404 0.0955 .
#> *--------
#> Parameter: sigma
#> ---
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.4196 0.3841 1.092 0.3545
#> x -2.2085 0.7490 -2.949 0.0601 .
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> *--------
#> n = 7 df = 4 res.df = 3
#> Deviance = 19.748 Null Dev. Red. = 20.82%
#> AIC = 27.748 elapsed = 0.02sec
# Obtaining the means
m1 <- fitted(mod2, model="mu", type="parameter")
m2 <- est_param_2(mod2, m=10000, fun="mean")
cbind(m1, m2) # Comparing the estimated means
#> m1 m2
#> 1 6.960302 6.849494
#> 2 1.388837 1.392931
#> 3 3.396289 3.392616
#> 4 2.932452 2.927701
#> 5 4.912643 4.919701
#> 6 1.213823 1.213344
#> 7 2.694837 2.708466
cor(m1, m2)
#> [1] 0.9999016
# Obtaining the standard deviations
s1 <- fitted(mod2, what="sigma") * fitted(mod2, what="mu")
s2 <- est_param_2(mod2, m=10000, fun="sd")
cbind(s1, s2) # Comparing the estimated standard deviations
#> s1 s2
#> 1 10.3652606 10.2748676
#> 2 0.3906638 0.3880437
#> 3 2.4083915 2.4229803
#> 4 1.7865251 1.8068101
#> 5 5.1027538 5.0345392
#> 6 0.2970436 0.2979131
#> 7 1.5044000 1.4955879
cor(s1, s2)
#> [1] 0.9999838
# To predict mean and standard deviation for new observations
new_data <- data.frame(x=c(0.26, 0.52, 0.77))
est_param_2(mod=mod2, fun="mean", m=10000, data=data2, newdata=new_data)
#> 1 2 3
#> 4.075023 2.351751 1.372751
est_param_2(mod=mod2, fun="sd" , m=10000, data=data2, newdata=new_data)
#> 1 2 3
#> 3.4902395 1.1234900 0.3825976