[This article was first published on Yet Another Blog in Statistical Computing » S+/R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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 ***
To leave a comment for the author, please follow the link and comment on their blog: Yet Another Blog in Statistical Computing » S+/R.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.