[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.
library(betareg)
library(sas7bdat)
df1 <- read.sas7bdat('lgd.sas7bdat')
df2 <- df1[df1$y < 1, ]
fml <- as.formula('y ~ x2 + x3 + x4 + x5 + x6 | x3 + x4 | x1 + x2')
### LATENT-CLASS BETA REGRESSION: AIC = -565 ###
mdl1 <- betamix(fml, data = df2, k = 2, FLXcontrol = list(iter.max = 500, minprior = 0.1))
print(mdl1)
#betamix(formula = fml, data = df2, k = 2, FLXcontrol = list(iter.max = 500,
# minprior = 0.1))
#
#Cluster sizes:
# 1 2
#157 959
summary(mdl1, which = 'concomitant')
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) -1.35153 0.41988 -3.2188 0.001287 **
#x1 2.92537 1.13046 2.5878 0.009660 **
#x2 2.82809 1.42139 1.9897 0.046628 *
summary(mdl1)
#$Comp.1$mean
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) -0.8963228 1.0385545 -0.8630 0.3881108
#x2 3.1769062 0.6582108 4.8266 1.389e-06 ***
#x3 -0.0520060 0.0743714 -0.6993 0.4843805
#x4 4.9642998 1.4204071 3.4950 0.0004741 ***
#x5 0.0021647 0.0022659 0.9554 0.3393987
#x6 0.0248573 0.0062982 3.9467 7.922e-05 ***
#
#$Comp.1$precision
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) -5.37817 1.44817 -3.7138 0.0002042 ***
#x3 0.45009 0.10094 4.4589 8.239e-06 ***
#x4 3.06969 1.41450 2.1702 0.0299948 *
#
#$Comp.2
#$Comp.2$mean
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) -1.8737088 0.3888454 -4.8186 1.445e-06 ***
#x2 -0.6318086 0.1892501 -3.3385 0.0008424 ***
#x3 0.1786425 0.0265428 6.7303 1.693e-11 ***
#x4 2.0646272 0.5256002 3.9281 8.561e-05 ***
#x5 -0.0064821 0.0014053 -4.6127 3.974e-06 ***
#x6 0.0018828 0.0022873 0.8231 0.4104318
#
#$Comp.2$precision
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) 1.092403 0.616974 1.7706 0.076630 .
#x3 0.017330 0.040024 0.4330 0.665029
#x4 2.138812 0.717702 2.9801 0.002882 **
### BETA REGRESSION TREE: AIC = -578 ###
mdl2 <- betatree(fml, data = df2, minsplit = 100)
print(mdl2)
#1) x2 <= 0.08584895; criterion = 1, statistic = 154.716
# 2)* weights = 121
#Terminal node model
#betaReg fit with coefficients:
# (Intercept) x2 x3 x4
# 3.307359 -2.854351 -0.262815 -2.414481
# x5 x6 (phi)_(Intercept) (phi)_x3
# -0.007555 0.030346 1.003767 -0.002907
# (phi)_x4
# 2.528602
#
#1) x2 > 0.08584895
# 3)* weights = 995
#Terminal node model
#betaReg fit with coefficients:
# (Intercept) x2 x3 x4
# -2.134931 -0.194830 0.168136 2.811077
# x5 x6 (phi)_(Intercept) (phi)_x3
# -0.002070 0.004677 -1.018102 0.151778
# (phi)_x4
# 2.142995
sctest(mdl2, node = 1)
# x1 x2
#statistic 113.4781 154.7165
#p.value 0.0000 0.0000
summary(mdl2)
#$`2`
#
#Coefficients (mean model with logit link):
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) 3.307359 1.091289 3.031 0.002440 **
#x2 -2.854351 3.644882 -0.783 0.433561
#x3 -0.262815 0.074716 -3.518 0.000436 ***
#x4 -2.414481 1.785447 -1.352 0.176276
#x5 -0.007555 0.002788 -2.710 0.006738 **
#x6 0.030346 0.006833 4.441 8.96e-06 ***
#
#Phi coefficients (precision model with log link):
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) 1.003767 1.353496 0.742 0.458
#x3 -0.002907 0.090816 -0.032 0.974
#x4 2.528602 2.344241 1.079 0.281
#$`3`
#
#Coefficients (mean model with logit link):
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) -2.134931 0.337784 -6.320 2.61e-10 ***
#x2 -0.194830 0.144062 -1.352 0.17625
#x3 0.168136 0.022521 7.466 8.28e-14 ***
#x4 2.811077 0.387788 7.249 4.20e-13 ***
#x5 -0.002070 0.001136 -1.822 0.06848 .
#x6 0.004677 0.001770 2.643 0.00822 **
#
#Phi coefficients (precision model with log link):
# Estimate Std. Error z value Pr(>|z|)
#(Intercept) -1.01810 0.46575 -2.186 0.028821 *
#x3 0.15178 0.03057 4.965 6.88e-07 ***
#x4 2.14300 0.56979 3.761 0.000169 ***
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.