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)
#> Loading required package: mgcv
#> This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
#>
#> Attaching package: 'mgcv'
#> The following object is masked from 'package:gamlss':
#>
#> lp
#>
#> Attaching package: 'gamlss2'
#> The following objects are masked from 'package:gamlss':
#>
#> GAIC, Rsq
f <- y ~ x | x
mod1 <- gamlss2(f, 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
f <- y ~ x | x
mod2 <- gamlss2(f, family=GA, data=data2)
#> GAMLSS-RS iteration 1: Global Deviance = 23.323 eps = 0.082109
#> GAMLSS-RS iteration 2: Global Deviance = 20.2365 eps = 0.132337
#> GAMLSS-RS iteration 3: Global Deviance = 19.751 eps = 0.023991
#> GAMLSS-RS iteration 4: Global Deviance = 19.748 eps = 0.000152
#> GAMLSS-RS iteration 5: Global Deviance = 19.748 eps = 0.000000
summary(mod2)
#> Call:
#> gamlss2(formula = f, data = data2, family = GA)
#> ---
#> Family: GA
#> Link function: mu = log, sigma = log
#> *--------
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> mu.(Intercept) 1.9621 0.6399 3.066 0.0547 .
#> mu.x -2.1376 0.8890 -2.404 0.0955 .
#> sigma.(Intercept) 0.4185 0.3838 1.090 0.3553
#> sigma.x -2.2057 0.7490 -2.945 0.0603 .
#> ---
#> 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.08sec
# 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.968609 6.857792
#> 2 1.388673 1.392771
#> 3 3.398358 3.394685
#> 4 2.933888 2.929135
#> 5 4.917113 4.924174
#> 6 1.213547 1.213068
#> 7 2.695972 2.709608
cor(m1, m2)
#> [1] 0.9999019
# 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.3662769 10.2760313
#> 2 0.3910305 0.3884071
#> 3 2.4095289 2.4241242
#> 4 1.7875057 1.8078031
#> 5 5.1041812 5.0359773
#> 6 0.2973433 0.2982147
#> 7 1.5052921 1.4964750
cor(s1, s2)
#> [1] 0.9999837
# 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.078112 2.352476 1.372576
est_param_2(mod=mod2, fun="sd" , m=10000, data=data2, newdata=new_data)
#> 1 2 3
#> 3.4915686 1.1242380 0.3829595
# Example with BCTo distribution ------------------------------------------
# Here E(y) = mu
n <- 7
x <- runif(n=n)
y <- rBCTo(n=n, mu=exp(-1 + 5 * x), sigma=exp(-1 + 2 * x))
data3 <- data.frame(y, x)
data3
#> y x
#> 1 34.053855 0.9407333
#> 2 1.099746 0.2202595
#> 3 41.438009 0.8612220
#> 4 148.279299 0.9834827
#> 5 8.928781 0.5866224
#> 6 2.032540 0.3853315
#> 7 17.896157 0.6776835
f <- y ~ x | x
mod3 <- gamlss2(f, family=BCTo, data=data3)
#> GAMLSS-RS iteration 1: Global Deviance = 46.0119 eps = 0.751454
#> GAMLSS-RS iteration 2: Global Deviance = 37.647 eps = 0.181797
#> GAMLSS-RS iteration 3: Global Deviance = 37.4287 eps = 0.005798
#> GAMLSS-RS iteration 4: Global Deviance = 37.2823 eps = 0.003913
#> GAMLSS-RS iteration 5: Global Deviance = 37.1781 eps = 0.002793
#> GAMLSS-RS iteration 6: Global Deviance = 37.0689 eps = 0.002938
#> GAMLSS-RS iteration 7: Global Deviance = 36.9562 eps = 0.003040
#> GAMLSS-RS iteration 8: Global Deviance = 36.8468 eps = 0.002959
#> GAMLSS-RS iteration 9: Global Deviance = 36.7366 eps = 0.002990
#> GAMLSS-RS iteration 10: Global Deviance = 36.6327 eps = 0.002828
#> GAMLSS-RS iteration 11: Global Deviance = 36.532 eps = 0.002748
#> GAMLSS-RS iteration 12: Global Deviance = 36.4347 eps = 0.002665
#> GAMLSS-RS iteration 13: Global Deviance = 36.3408 eps = 0.002577
#> GAMLSS-RS iteration 14: Global Deviance = 36.2491 eps = 0.002520
#> GAMLSS-RS iteration 15: Global Deviance = 36.1611 eps = 0.002427
#> GAMLSS-RS iteration 16: Global Deviance = 36.0748 eps = 0.002386
#> GAMLSS-RS iteration 17: Global Deviance = 35.9899 eps = 0.002354
#> GAMLSS-RS iteration 18: Global Deviance = 35.9078 eps = 0.002282
#> GAMLSS-RS iteration 19: Global Deviance = 35.8263 eps = 0.002269
#> GAMLSS-RS iteration 20: Global Deviance = 35.7465 eps = 0.002226
summary(mod3)
#> Call:
#> gamlss2(formula = f, data = data3, family = BCTo)
#> ---
#> Family: BCTo
#> Link function: mu = log, sigma = log, nu = identity, tau = log
#> *--------
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> mu.(Intercept) -1.095 1.003 -1.092 0.472
#> mu.x 5.308 3.735 1.421 0.390
#> sigma.(Intercept) -4.150 8.964 -0.463 0.724
#> sigma.x 5.293 21.226 0.249 0.844
#> nu.(Intercept) 2.896 5.112 0.567 0.672
#> tau.(Intercept) 15.896 60.012 0.265 0.835
#> *--------
#> n = 7 df = 6 res.df = 1
#> Deviance = 35.7465 Null Dev. Red. = 43.24%
#> AIC = 47.7465 elapsed = 0.74sec
# Obtaining the means
m1 <- fitted(mod3, model="mu", type="parameter")
m2 <- est_param_2(mod3, m=10000, fun="mean")
cbind(m1, m2) # Comparing the estimated means
#> m1 m2
#> 1 49.302617 83.042891
#> 2 1.076669 1.074539
#> 3 32.328534 47.909537
#> 4 61.860459 111.952411
#> 5 7.526572 7.799242
#> 6 2.585801 2.550140
#> 7 12.204086 13.875374
cor(m1, m2)
#> [1] 0.9959884
# To predict mean and for new observations
new_data <- data.frame(x=c(0.26, 0.52, 0.77))
est_param_2(mod=mod3, fun="mean", m=10000, data=data3, newdata=new_data)
#> 1 2 3
#> 1.325256 5.244328 25.414443