## Estimating a Beta Regression with The Variable Dispersion in R

pkgs <- c('sas7bdat', 'betareg', 'lmtest') lapply(pkgs, require, character.only = T) df1 <- read.sas7bdat("lgd.sas7bdat") df2 <- df1[which(df1$y < 1), ] xvar <- paste("x", 1:7, sep = '', collapse = " + ") fml1 <- as.formula(paste("y ~ ", xvar)) fml2 <- as.formula(paste("y ~ ", xvar, "|", xvar)) # FIT A BETA MODEL WITH THE FIXED PHI beta1 <- betareg(fml1, data = df2) summary(beta1) # Coefficients (mean model with logit link): # Estimate Std. Error z value Pr(>|z|) # (Intercept) -1.500242 0.329670 -4.551 5.35e-06 *** # x1 0.007516 0.026020 0.289 0.772680 # x2 0.429756 0.135899 3.162 0.001565 ** # x3 0.099202 0.022285 4.452 8.53e-06 *** # x4 2.465055 0.415657 5.931 3.02e-09 *** # x5 -0.003687 0.001070 -3.446 0.000568 *** # x6 0.007181 0.001821 3.943 8.06e-05 *** # x7 0.128796 0.186715 0.690 0.490319 # # Phi coefficients (precision model with identity link): # Estimate Std. Error z value Pr(>|z|) # (phi) 3.6868 0.1421 25.95 <2e-16 *** # FIT A BETA MODEL WITH THE VARIABLE PHI beta2 <- betareg(fml2, data = df2) summary(beta2) # Coefficients (mean model with logit link): # Estimate Std. Error z value Pr(>|z|) # (Intercept) -1.996661 0.336445 -5.935 2.95e-09 *** # x1 0.007033 0.026809 0.262 0.793072 # x2 0.371098 0.135186 2.745 0.006049 ** # x3 0.133356 0.022624 5.894 3.76e-09 *** # x4 2.951245 0.401493 7.351 1.97e-13 *** # x5 -0.003475 0.001119 -3.105 0.001902 ** # x6 0.006528 0.001884 3.466 0.000529 *** # x7 0.100274 0.190915 0.525 0.599424 # # Phi coefficients (precision model with log link): # Estimate Std. Error z value Pr(>|z|) # (Intercept) -0.454399 0.452302 -1.005 0.315072 # x1 0.009119 0.035659 0.256 0.798150 # x2 0.611049 0.188225 3.246 0.001169 ** # x3 0.092073 0.030678 3.001 0.002689 ** # x4 2.248399 0.579440 3.880 0.000104 *** # x5 -0.002188 0.001455 -1.504 0.132704 # x6 -0.000317 0.002519 -0.126 0.899847 # x7 -0.166226 0.256199 -0.649 0.516457 # LIKELIHOOD RATIO TEST TO COMPARE BOTH BETA MODELS lrtest(beta1, beta2) # Likelihood ratio test # # Model 1: y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 # Model 2: y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 | x1 + x2 + x3 + x4 + x5 + x6 + x7 # #Df LogLik Df Chisq Pr(>Chisq) # 1 9 231.55 # 2 16 257.24 7 51.38 7.735e-09 ***

Advertisements