R/est_param.R
est_param.RdThis function estimates parameters like mean, standard deviation for a gamlss object. This function is very useful when the mean or standard deviation do not match with mu, sigma or some function of mu and/or sigma.
est_param(mod, fun = "mean", m = 100, data, newdata = NULL)a gamlss object.
the estimated parameter, by default is mean but can be change for sd or var.
a number to indicate the number of observations to simulate, by default its values is 100.
the original data frame used to fit mod.
a data frame containing new values for the explanatory variables.
est_param 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 -5.962582 0.080750138
#> 2 14.118155 0.834333037
#> 3 -7.442000 0.600760886
#> 4 -5.164947 0.157208442
#> 5 -7.290830 0.007399441
#> 6 -10.870123 0.466393497
#> 7 14.465420 0.497777389
library(gamlss)
#> Loading required package: splines
#> Loading required package: gamlss.data
#>
#> Attaching package: 'gamlss.data'
#> The following object is masked from 'package:datasets':
#>
#> sleep
#> Loading required package: gamlss.dist
#> Loading required package: nlme
#> Loading required package: parallel
#> ********** GAMLSS Version 5.4-20 **********
#> For more on GAMLSS look at https://www.gamlss.com/
#> Type gamlssNews() to see new features/changes/bug fixes.
mod1 <- gamlss(y ~ x, sigma.fo=~x, family=NO, data=data1)
#> GAMLSS-RS iteration 1: Global Deviance = 43.4364
#> GAMLSS-RS iteration 2: Global Deviance = 35.2613
#> GAMLSS-RS iteration 3: Global Deviance = 34.1276
#> GAMLSS-RS iteration 4: Global Deviance = 34.0948
#> GAMLSS-RS iteration 5: Global Deviance = 34.0943
# Obtaining the means
m1 <- fitted(mod1, what="mu")
m2 <- est_param(mod1, m=10000, fun="mean")
cbind(m1, m2) # Comparing the estimated means
#> m1 m2
#> 1 -6.09961525 -6.0990137
#> 2 5.96700504 9.0397923
#> 3 2.22696999 2.6019221
#> 4 -4.87533943 -4.8755479
#> 5 -7.27413102 -7.2749168
#> 6 0.07543462 0.1363465
#> 7 0.57796394 0.4918223
cor(m1, m2)
#> [1] 0.9899216
# Obtaining the standard deviations
s1 <- fitted(mod1, what="sigma")
s2 <- est_param(mod1, m=10000, fun="sd")
cbind(s1, s2) # Comparing the estimated standard deviations
#> s1 s2
#> 1 0.13810096 0.13986106
#> 2 276.18534671 274.26679377
#> 3 26.18534313 26.31902698
#> 4 0.29861774 0.29864206
#> 5 0.06590076 0.06603477
#> 6 6.75264134 6.75616393
#> 7 9.26719010 9.29111948
cor(s1, s2)
#> [1] 0.9999994
# Example for GAmma regression --------------------------
# Here E(y) = mu and Var(y)= sigma^2 mu^2
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.04439329 0.93620147
#> 2 24.07903452 0.49056477
#> 3 0.77665476 0.08767044
#> 4 0.19226774 0.90254377
#> 5 1.66785029 0.27514128
#> 6 0.41760508 0.18782978
#> 7 380.30685933 0.98977368
mod2 <- gamlss(y ~ x, sigma.fo=~x, family=GA, data=data2)
#> GAMLSS-RS iteration 1: Global Deviance = 34.2295
#> GAMLSS-RS iteration 2: Global Deviance = 33.7537
#> GAMLSS-RS iteration 3: Global Deviance = 33.7535
summary(mod2)
#> Warning: summary: vcov has failed, option qr is used instead
#> ******************************************************************
#> Family: c("GA", "Gamma")
#>
#> Call: gamlss(formula = y ~ x, sigma.formula = ~x, family = GA, data = data2)
#>
#>
#> Fitting method: RS()
#>
#> ------------------------------------------------------------------
#> Mu link function: log
#> Mu Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -1.1351 0.4518 -2.512 0.05367 .
#> x 6.9475 1.5450 4.497 0.00642 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ------------------------------------------------------------------
#> Sigma link function: log
#> Sigma Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.9437 0.4610 -2.047 0.096 .
#> x 1.9606 0.6402 3.063 0.028 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> ------------------------------------------------------------------
#> No. of observations in the fit: 7
#> Degrees of Freedom for the fit: 4
#> Residual Deg. of Freedom: 3
#> at cycle: 3
#>
#> Global Deviance: 33.75354
#> AIC: 41.75354
#> SBC: 41.53718
#> ******************************************************************
# Obtaining the means
m1 <- fitted(mod2, what="mu")
m2 <- est_param(mod2, m=10000, fun="mean")
cbind(m1, m2) # Comparing the estimated means
#> m1 m2
#> 1 214.6680721 220.8212299
#> 2 9.7089855 9.6087059
#> 3 0.5909412 0.5875801
#> 4 169.9080886 175.1053139
#> 5 2.1736475 2.1644163
#> 6 1.1850842 1.1858251
#> 7 311.4639353 294.1534341
cor(m1, m2)
#> [1] 0.9984612
# Obtaining the standard deviations
s1 <- fitted(mod2, what="sigma") * fitted(mod2, what="mu")
s2 <- est_param(mod2, m=10000, fun="sd")
cbind(s1, s2) # Comparing the estimated standard deviations
#> s1 s2
#> 1 523.6485741 524.2365047
#> 2 9.8856690 10.1036040
#> 3 0.2731049 0.2753241
#> 4 387.9970826 371.6552429
#> 5 1.4507668 1.4474450
#> 6 0.6665239 0.6701837
#> 7 843.9071060 813.6235151
cor(s1, s2)
#> [1] 0.9997663
# To predict mean and standard deviation for new observations
new_data <- data.frame(x=c(0.26, 0.52, 0.77))
est_param(mod=mod2, fun="mean", m=10000, data=data2, newdata=new_data)
#> [1] 1.955015 11.809437 68.535349
est_param(mod=mod2, fun="sd" , m=10000, data=data2, newdata=new_data)
#> [1] 1.269137 12.903712 117.566921