Two-Part Model with Marginal Effects.
knitr::opts_chunk$set(echo = TRUE)
library(twopartm)
1. Continuous health cost data example: MEPS
Data
##data about health expenditures, i.e., non-negative continuous response
data(meps,package = "twopartm")
##Data information
?meps
Fit two-part model:
Fit two-part model with the same regressors in both parts, with logistic regression model for the first part, and glm with Gamma family with log link for the second-part model
tpmodel = tpm(exp_tot~female+age, data = meps,link_part1 = "logit",family_part2 = Gamma(link = "log"))
tpmodel
fit two-part model with different regressors in both parts
tpmodel = tpm(formula_part1 = exp_tot~female+age, formula_part2 = exp_tot~female+age+ed_colplus,data = meps,link_part1 = "logit",family_part2 = Gamma(link = "log"))
tpmodel
fit two-part model with transformed regressors and randomly assigned weights
set.seed(100)
meps$weights = sample(1:30,nrow(meps),replace = T)
tpmodel = tpm(formula_part1 = exp_tot~female+age, formula_part2 = exp_tot~female+I(age^2)+ed_colplus,data = meps,link_part1 = "logit",family_part2 = Gamma(link = "log"),weights = meps$weights)
tpmodel
Model object
##fit two-part model with the same regressors in both parts
tpmodel = tpm(exp_tot~female+age, data = meps,link_part1 = "logit",family_part2 = Gamma(link = "log"))
tpmodel
Get the formula specified for the first-part model
tpmodel@formula_part1
Get the formula specified for the second-part model
tpmodel@formula_part2
Get the log-likelihood for the fitted two-part model
tpmodel@loglik
Get the the fitted glm model for the first part
tpmodel@model_part1
Get the the fitted glm model for the second part
tpmodel@model_part2
Methods about two-part model object
##information about fitted two-part model
print(tpmodel)
##summary information
summary(tpmodel)
##estimated coefficients for both parts
coef(tpmodel)
##estimated coefficients for the first-part model
coef(tpmodel,model = "model1")
##response residues from the full two-part model
res = residuals(tpmodel)
##response residues from the first-part model
res1 = residuals(tpmodel,model = "model1")
##deviance residues from the second-part model
res2 = residuals(tpmodel,model = "model2",type = "deviance")
##log-likehood
logLik(tpmodel)
##plots for two-part model
plot(tpmodel)
Prediction
Get prediction results with standard errors for the first 10 observations in the dataset
predict(tpmodel,newdata = meps[1:10,],se.fit = T)
Average Marginal Effect (AME)
Fit two-part model with different regressors in both parts
##fit two-part model with different regressors in both parts
tpmodel = tpm(formula_part1 = exp_tot~female+age, formula_part2 = exp_tot~female+age+ed_colplus,data = meps,link_part1 = "logit",family_part2 = Gamma(link = "log"))
tpmodel
AMEs for all variables with standard errors and CIs
AME(tpmodel)
AMEs for variable "female"
AME(tpmodel,term = "female")
AMEs for variables "female" and "age" with standard errors by bootstrap methods, CIs by bootstrap quantiles at level 0.9
AME(tpmodel,term = c("female","age"),se.method = "bootstrap", CI.boots = T, level = 0.9)
AMEs for all variables with standard errors and CIs at age 20,40,60,80 respectively
AME(tpmodel,at = list(age = c(20,40,60,80)))
AMEs for all variables with standard errors and CIs at age 20,40,60,80, and education level is more than college
AME(tpmodel,at = list(age = c(20,40,60,80),ed_colplus = 1),term = "female",se.method = "bootstrap")
Predictive Margins with Ratios
Predictive margins and corresponding ratios for all variables with standard errors and CIs. For factor or logical variables, predictive margins at all the levels are calculated, and for numeric (and integer) variables, predictive margins at the mean values among observations are calculated.
margin(tpmodel)
Predictive margins and corresponding ratios for variable "age" at 20,40,60,80, with standard errors and CIs.
margin(tpmodel,term = "age",value = list(age = c(20,40,60,80)))
Predictive margins and corresponding ratios for female, age at 20,40,60,80, and more than college education level, respectively
margin(tpmodel,value = list(female = 1,age = c(50,70),ed_colplus = 1))
Predictive margins and corresponding ratios for variable "ed_colplus" with standard errors by bootstrap methods, and CIs by bootstrap quantiles at level 0.99
margin(tpmodel,term = "ed_colplus",se.method = "bootstrap",CI.boots = T, level = 0.99)
Predictive margins and corresponding ratios for all variables with standard errors and CIs calculated on the first 500 observations
margin(tpmodel,newdata = meps[1:500,])
2. Count data example: bioChemistry
Data
data("bioChemists",package = "twopartm")
##Data information
?bioChemists
Fit two-part model:
Fit two-part model with the same regressors in both parts, with logistic regression model for the first part, and Poisson regression model with default log link for the second-part model
tpmodel = tpm(art ~ .,data = bioChemists,link_part1 = "logit",family_part2 = poisson)
tpmodel
summary(tpmodel)
##estimated coefficients for both parts
coef(tpmodel)
##log-likehood
logLik(tpmodel)
##plots for two-part model
plot(tpmodel)
Average Marginal Effect (AME)
AMEs for all variables with standard errors and CIs
AME(tpmodel)
AMEs for variables "fem" and "kid5" with standard errors
AME(tpmodel,term = c("fem","kid5"))
AMEs for all variables if all are women
AME(tpmodel,at = list(fem = "Women"))
AMEs for variable "ment" when all are women and the numbers of children aged 5 or younger are 0,1,3,5, with standard errors by bootstrap methods, and CIs by bootstrap quantiles
AME(tpmodel,term = "ment",at = list(fem = "Women",kid5 = c(0,1,3,5)),se.method = "bootstrap",CI.boots = T)
Predictive Margins with Ratios
Predictive margins and corresponding ratios for all variables with standard errors and CIs.
margin(tpmodel)
Predictive margins and corresponding ratios for variable "kid5" at 1,2,3,5, with standard errors by bootstrap methods, and CIs by bootstrap quantiles
margin(tpmodel,term = "kid5",value = list(kid5 = c(1,2,3,5)),se.method = "bootstrap",CI.boots = T)
Predictive margins and corresponding ratios for variable "ment" at 6,7,8, without standard errors and CIs
margin(tpmodel,term = "ment",value = list(ment = c(6,7,8)),se = F)
Predictive margins and corresponding ratios for women and all the levels of variable "mar",with standard errors by bootstrap methods, and normal-based CIs
margin(tpmodel,term = c("fem","mar"),value = list(fem = "Women"),se.method = "bootstrap")
Predictive margins and corresponding ratios for all the levels of variable "mar", and for variable "phd" at 2.5,3,3.2, calculated on the first 500 observations, with standard errors and CIs
margin(tpmodel,newdata = bioChemists[1:500,],term = c("phd","mar"),value = list(phd = c(2.5,3,3.2)))