Modelos Tobit y en dos partes

15/4/2021

DATA del Medical Expenditure Panel Survey (MEPS)

#Lectura de datos
Hexpend = read.csv('HealthExpendEvent.csv',  quote = "",header=TRUE)
#  Checar nombres, dimnesión y primeras 8 observaciones ;
names(Hexpend)
dim(Hexpend)
Hexpend[1:8,]
attach(Hexpend)
#  Recodificación y exploración de variables ;
require(ggplot2)
ggplot(Hexpend,aes(x = RACE))+geom_bar()+coord_flip()
ASIAN=1*(RACE1==1)
BLACK=1*(RACE1==2)
NATIVE=1*(RACE1==3)
WHITE=1-ASIAN-BLACK-NATIVE   # combinar other y white

ggplot(Hexpend,aes(x = REGION))+geom_bar()+coord_flip()
NORTHEAST=1*(REGION1==1)
MIDWEST=1*(REGION1==2)
SOUTH=1*(REGION1==3)
WEST=1-NORTHEAST-MIDWEST-SOUTH

ggplot(Hexpend,aes(x = EDUC))+geom_bar()+coord_flip()
HIGHSCHOOL=1*(EDUC1==1)
COLLEGE=1*(EDUC1==2)
LOWERHS=1-HIGHSCHOOL-COLLEGE

ggplot(Hexpend,aes(x = PHSTAT))+geom_bar()+coord_flip()
POOR=1*(PHSTAT1==4)
FAIR=1*(PHSTAT1==3)
GOOD=1*(PHSTAT1==2)
VGOOD=1*(PHSTAT1==1)
EXCELLNT=1-POOR-FAIR-GOOD-VGOOD

ggplot(Hexpend,aes(x = INCOME))+geom_bar()+coord_flip()
NPOOR=1*(INCOME1==1)
LINCOME=1*(INCOME1==2)
MINCOME=1*(INCOME1==3)
HINCOME=1*(INCOME1==4)
POORNEG=1-NPOOR-LINCOME-MINCOME-HINCOME
ggplot(Hexpend,aes(x = INCOME,fill=PHSTAT))+geom_bar()+coord_flip()
GENDER1 <- as_factor(Hexpend$GENDER)
levels(GENDER1) <- c("MALE","FEMALE")
ggplot(Hexpend,aes(x = INCOME,fill=GENDER1))+geom_bar()+coord_flip()
ggplot(Hexpend,aes(x = GENDER1,fill=PHSTAT))+geom_bar()+coord_flip()
ggplot(Hexpend,aes(x = EDUC,fill=GENDER1))+geom_bar()+coord_flip()
ggplot(Hexpend,aes(x = INCOME,fill=PHSTAT))+geom_bar()+coord_flip()

COUNT=COUNTIP
EXP=EXPENDIP
POSEXP=1*(EXP>0)
LNEXP=ifelse(EXP>0,log(EXP),0)
LOWERBOUND=ifelse(EXP>0, LNEXP,'')

RaceFactor=cbind(ASIAN,BLACK,NATIVE,WHITE)
RegionFactor=cbind(NORTHEAST,MIDWEST,SOUTH,WEST)
EducFactor=cbind(HIGHSCHOOL,COLLEGE,LOWERHS)
PhstatFactor=cbind(POOR,FAIR,GOOD,VGOOD,EXCELLNT)
IncomeFactor=cbind(NPOOR,LINCOME,MINCOME,HINCOME,POORNEG)
PersonLevel=cbind(Hexpend[,-4],RaceFactor,RegionFactor,EducFactor,
                  PhstatFactor,IncomeFactor,COUNT,EXP,POSEXP,LNEXP,
                  LOWERBOUND)
names(PersonLevel)
dim(PersonLevel)

# Gastos positivos ;
PersonLevelPOSEXP=subset(PersonLevel,POSEXP==1)
names(PersonLevelPOSEXP)
dim(PersonLevelPOSEXP)
PersonLevelPOSEXP[1:8,]

Estadísticas de resumen

require(Hmisc) #install.packages("Hmisc")  "Harrel-verse"  
attach(PersonLevel)
count=as.data.frame(table(GENDER)) #  Por Género;
colnames(count)=c("GENDER","TotalCount")
count
meanPOS=summarize(POSEXP, GENDER, mean)
meanPOS
meanLNEXP=summarize(LNEXP,GENDER,mean)
meanLNEXP                                  
detach(PersonLevel)
attach(PersonLevelPOSEXP)
POScount=as.data.frame(table(GENDER))
POScount
POSmeanLNEXP=summarize(PersonLevelPOSEXP$LNEXP,
                       PersonLevelPOSEXP$GENDER,mean)
POSmeanLNEXP

Table1Gender1=cbind(count,meanPOS,meanLNEXP,POScount,POSmeanLNEXP)
Table1Gender1$GENDER=NULL
Table1Gender1$GENDER=NULL
Table1Gender1$GENDER=NULL
Table1Gender1$GENDER=NULL
attach(Table1Gender1)
Percent=Table1Gender1$TotalCount/2000 
#se necesita el número total de observaciones;
Variable=Table1Gender1$GENDER
Label=c("MALE","FEMALE")
Table1Gender1$GENDER=NULL
Table1Gender1$TotalCount=NULL
Table1Gender1$Freq=NULL
Table1Gender=cbind(Percent,Table1Gender1,Variable,Label)
colnames(Table1Gender)=c("Percent","MeanPosExp","MeanLnExp",
                         "PosMean","Variable","Label")
Table1Gender

# Automatización del proceso ;
SumStat<-function(x,y)
{
  attach(PersonLevel)
  count=as.data.frame(table(x))
  colnames(count)=c("Vari","TotalCount")
  meanPOS=summarize(PersonLevel$POSEXP,x,mean)
  colnames(meanPOS)=c("Vari","meanPOS")
  meanLNEXP=summarize(LNEXP,x,mean)
  colnames(meanLNEXP)=c("Vari","LNEXP")
  POScount=as.data.frame(table(y))
  colnames(POScount)=c("Vari","Freq")
  POSmeanLNEXP=summarize(PersonLevelPOSEXP$LNEXP,y,mean)
  colnames(POSmeanLNEXP)=c("Vari","PersonLevelPOSEXP$LNEXP")

# Automatización (cont.)
  Table1y1=cbind(count,meanPOS,meanLNEXP,POScount,POSmeanLNEXP)
  Table1y1$Vari=NULL
  Table1y1$Vari=NULL
  Table1y1$Vari=NULL
  Table1y1$Vari=NULL
  Percent=Table1y1$TotalCount/2000  
  Variable=Table1y1$Vari                        
  Table1y1$Vari=NULL
  Table1y1$TotalCount=NULL
  Table1y1$Freq=NULL
  Table1y=cbind(Percent,Table1y1,Variable)
  colnames(Table161y)=c("Percent","MeanPosExp","MeanLnExp","PosMean",
                        "Variable")
  Table1y
}

# Aplicacamos SumStat()
TabGender=SumStat(GENDER,PersonLevelPOSEXP$GENDER)
Label=c("MALE","FEMALE")
Table1Gender=cbind(TabGender,Label)
Table1Gender
ggplot(Hexpend,aes(x = AGE,y=LNEXP,colour=GENDER1))+geom_point()

TabRace=SumStat(RACE1,PersonLevelPOSEXP$RACE1)
Label=c("OTHERS","ASIAN","BLACK","NATIVE","WHITE")
Table1Race=cbind(TabRace,Label)
Table1Race
ggplot(Hexpend,aes(x = AGE,y=LNEXP,shape=RACE))+geom_point()

TabRegion=SumStat(REGION1,PersonLevelPOSEXP$REGION1)
Label=c("WEST","NORTHEAST","MIDWEST","SOUTH")
Table1Region=cbind(TabRegion,Label)
Table1Region
ggplot(Hexpend,aes(x = AGE,y=LNEXP,colour=REGION))+geom_point()

TabEduc=SumStat(EDUC1,PersonLevelPOSEXP$EDUC1)
Label=c("LOWERHS","HIGHSCHOOL","COLLEGE")
Table1Educ=cbind(TabEduc,Label)
Table1Educ
ggplot(Hexpend,aes(x = AGE,y=LNEXP,colour=EDUC))+geom_point()

TabPhstat=SumStat(PHSTAT1,PersonLevelPOSEXP$PHSTAT1)
Label=c("EXCELLENT","VGOOD","GOOD","FAIR","POOR")
Table1Phstat=cbind(TabPhstat,Label)
Table1Phstat
ggplot(Hexpend,aes(x = AGE,y=LNEXP,shape=PHSTAT))+geom_point()

TabMnh=SumStat(MNHPOOR,PersonLevelPOSEXP$MNHPOOR)
Label=c("GOOD","POOR")
Table1Mnh=cbind(TabMnh,Label)
Table1Mnh
ggplot(Hexpend,aes(x = MNHPOOR,y=LNEXP,colour=GENDER1))+geom_point()

TabAnylimit=SumStat(ANYLIMIT,PersonLevelPOSEXP$ANYLIMIT)
Label=c("NO","YES")
Table1Anylimit=cbind(TabAnylimit,Label)
Table1Anylimit
ggplot(Hexpend,aes(x = ANYLIMIT,y=LNEXP,colour=GENDER1))+geom_point()

TabIncome=SumStat(INCOME1,PersonLevelPOSEXP$INCOME1)
Label=c("POOR","NPOOR","LINCOME","MINCOME","HINCOME")
Table1Income=cbind(TabIncome,Label)
Table1Income
ggplot(Hexpend,aes(x = AGE,y=LNEXP,shape=INCOME))+geom_point()

TabInsure=SumStat(INSURE,PersonLevelPOSEXP$INSURE)
Label=c("UNINSURED","INSURED")
Table1Insure=cbind(TabInsure,Label)
Table1Insure
ggplot(Hexpend,aes(x = INSURE,y=LNEXP,colour=GENDER1))+geom_point()

# Crear tabla completa
Table1n=rbind(Table1Gender,Table1Race,Table1Region,Table1Educ,
              Table1Phstat,Table1Mnh,Table1Anylimit,Table1Income,
              Table1Insure)
Category=c("DEMOGRAPHY"," ","ETHNICITY"," "," "," "," ","REGION"," ",
           " "," ","EDUCATION"," "," ","PHSTAT"," "," "," "," ", 
           "MNHEALTH"," ","ANYLIMIT"," ", "INCOME"," "," "," "," ", 
           "INSURED", " ")
Category
Table161=cbind(Table161n,Category)
require(pander)
pander(Table1)
write.csv(Table1,file="Table1")

Modelación

require(tidymodels)
require(dotwhisker)
#  Tabla Regresión OLS ;
lm_mod <- linear_reg() %>% 
  set_engine("lm")
lm_fit <- lm_mod %>% 
    fit(LNEXP ~ AGE+GENDER+ASIAN+BLACK+NATIVE+NORTHEAST+MIDWEST+SOUTH
       +COLLEGE+HIGHSCHOOL+POOR+FAIR+GOOD+VGOOD+MNHPOOR+ANYLIMIT
       +HINCOME+MINCOME+LINCOME+NPOOR+INSURE,data=PersonLevel)
lm_fit
broom::tidy(lm_fit)
broom::tidy(lm_fit) %>% 
  dwplot(dot_args = list(size = 2, color = "black"),
         whisker_args = list(color = "black"),
         vline=geom_vline(xintercept=0,colour='grey50',linetype=2))

#  Regresión Tobit: definición como regresión censurada;
require(survival)
surv_mod<-surv_reg(dist='gaussian') %>%
  set_engine("survival") 
surv_fit <- surv_mod %>% 
    fit(Surv(LNEXP,LNEXP>0,type="left") ~ AGE+GENDER+ASIAN+BLACK
       +NATIVE+NORTHEAST+MIDWEST+SOUTH
       +COLLEGE+HIGHSCHOOL+POOR+FAIR+GOOD+VGOOD+MNHPOOR+ANYLIMIT
       +HINCOME+MINCOME+LINCOME+NPOOR+INSURE,data=PersonLevel)
broom::tidy(surv_fit)
broom::tidy(surv_fit) %>% 
  dwplot(dot_args = list(size = 2, color = "black"),
         whisker_args = list(color = "black"),
         vline=geom_vline(xintercept=0,colour='grey50',linetype=2))

require(VGAM) #  Regresión Tobit empleando funciones específicas ;
summary(m <- vglm(LNEXP ~ AGE+GENDER+ASIAN+BLACK+NATIVE+NORTHEAST
                  +MIDWEST+SOUTH+COLLEGE+HIGHSCHOOL+POOR+FAIR+GOOD
                  +VGOOD+MNHPOOR+ANYLIMIT+HINCOME+MINCOME+LINCOME
                  +NPOOR+INSURE, tobit(Lower = 0, Upper=Inf), 
                  data = PersonLevel))
require(AER)
summary(m1<-tobit(LNEXP ~ AGE+GENDER+ASIAN+BLACK+NATIVE+NORTHEAST
                  +MIDWEST+SOUTH+COLLEGE+HIGHSCHOOL+POOR+FAIR+GOOD
                  +VGOOD+MNHPOOR+ANYLIMIT+HINCOME+MINCOME+LINCOME
                  +NPOOR+INSURE, left=0,right=Inf,dist="gaussian",
                  data = PersonLevel))
require(mgcv);require(cenGAM)
summary(m2 <- gam(LNEXP ~ AGE+GENDER+ASIAN+BLACK+NATIVE+NORTHEAST
                  +MIDWEST+SOUTH+COLLEGE+HIGHSCHOOL+POOR+FAIR+GOOD
                  +VGOOD+MNHPOOR+ANYLIMIT+HINCOME+MINCOME+LINCOME
                  +NPOOR+INSURE, family = tobit1(left.threshold=0)
                  ,data = PersonLevel))
m2$family$getTheta(TRUE) # sigma=exp(theta)

# Algoritmo de HECKMAN en dos etapas ; 
# Etapa 1 - Modelo Probit
glm_mod<-logistic_reg() %>% 
    set_engine("glm",family=binomial(link='probit')) %>% 
      set_mode("classification") 
glm_fit <- glm_mod %>% 
  fit(as.factor(POSEXP) ~ AGE+GENDER+ASIAN+BLACK+NATIVE+NORTHEAST
      +MIDWEST+SOUTH+COLLEGE+HIGHSCHOOL+POOR+FAIR+GOOD+VGOOD+MNHPOOR
      +ANYLIMIT+HINCOME+MINCOME+LINCOME+NPOOR+INSURE,data=PersonLevel)
broom::tidy(glm_fit)
broom::tidy(glm_fit) %>% 
  dwplot(dot_args = list(size = 2, color = "black"),
         whisker_args = list(color = "black"),
         vline=geom_vline(xintercept=0,colour="grey50",linetype=2))
# Relación Inversa de Mill
invMillRatio=dnorm(model1S$fitted.values)/pnorm(model1S$fitted.values)
ordernum=seq(1,2046)
PLnew=cbind(PersonLevel,invMillRatio,ordernum)

# Algoritmo de HECKMAN
# Etapa 2 
lm_mod2 <- linear_reg() %>% 
    set_engine("lm")
lm_fit2 <- lm_mod2 %>% 
    fit(LNEXP ~ AGE+GENDER+ASIAN+BLACK+NATIVE+NORTHEAST+MIDWEST+SOUTH
      +COLLEGE+HIGHSCHOOL+POOR+FAIR+GOOD+VGOOD+MNHPOOR+ANYLIMIT
      +HINCOME+MINCOME+LINCOME+NPOOR+INSURE+invMillRatio,data=PLnew)
lm_fit2
broom::tidy(lm_fit2)
broom::tidy(lm_fit2) %>% 
  dwplot(dot_args = list(size = 2, color = "black"),
         whisker_args = list(color = "black"),
         vline=geom_vline(xintercept=0,colour="grey50",linetype=2))

# Modelo en dos partes 
# Etapa 1;
glm_mod2<-logistic_reg() %>% 
    set_engine("glm",family=binomial(link='probit')) %>% 
      set_mode("classification") 
glm_fit2 <- glm_mod2 %>% 
    fit(as.factor(POSEXP) ~ AGE+GENDER+ASIAN+BLACK+NATIVE+NORTHEAST
        +MIDWEST+SOUTH+COLLEGE+HIGHSCHOOL+POOR+FAIR+GOOD+VGOOD
        +MNHPOOR+ANYLIMIT+HINCOME+MINCOME+LINCOME+NPOOR+INSURE,
        data=PersonLevel)
broom::tidy(glm_fit2)
broom::tidy(glm_fit2) %>% 
  dwplot(dot_args = list(size = 2, color = "black"),
         whisker_args = list(color = "black"),
         vline=geom_vline(xintercept=0,colour="grey50",linetype=2))

# Modelo en dos partes  
# Etapa 2;
glm_mod3<-logistic_reg() %>% 
    set_engine("glm",family=gaussian(link='identity')) %>% 
      set_mode("regression") 
glm_fit3 <- glm_mod3 %>% 
  fit(LNEXP ~ AGE+GENDER+ASIAN+BLACK+NATIVE+NORTHEAST+MIDWEST+SOUTH
      +COLLEGE+HIGHSCHOOL+POOR+FAIR+GOOD+VGOOD+MNHPOOR+ANYLIMIT
      +HINCOME+MINCOME+LINCOME+NPOOR+INSURE,data=PersonLevel)
broom::tidy(glm_fit3)
broom::tidy(glm_fit3) %>% 
  dwplot(dot_args = list(size = 2, color = "black"),
         whisker_args = list(color = "black"),
         vline=geom_vline(xintercept=0,colour="grey50",linetype=2))

# Modelo en dos partes 
# Etapa 1 - Variables significativas;
glm_mod4<-logistic_reg() %>% 
    set_engine("glm",family=binomial(link='probit')) %>% 
      set_mode("classification") 
glm_fit4 <- glm_mod4 %>% 
    fit(LNEXP ~ AGE+GENDER+ASIAN+BLACK+NATIVE+NORTHEAST+MIDWEST+SOUTH
      +COLLEGE+HIGHSCHOOL+POOR+FAIR+GOOD+VGOOD+MNHPOOR+ANYLIMIT
      +HINCOME+MINCOME+LINCOME+NPOOR+INSURE,data=PersonLevel)
broom::tidy(glm_fit4)
broom::tidy(glm_fit4) %>% 
  dwplot(dot_args = list(size = 2, color = "black"),
         whisker_args = list(color = "black"),
         vline=geom_vline(xintercept=0,colour="grey50",linetype=2))

# Modelo en dos partes 
# Etapa 2 - Variables significativas;
glm_mod5<-logistic_reg() %>% 
    set_engine("glm",family=gaussian) %>% 
      set_mode("regression") 
glm_fit5 <- glm_mod5 %>% 
  fit(LNEXP ~ AGE+GENDER+ASIAN+BLACK+NATIVE+NORTHEAST+MIDWEST+SOUTH
      +COLLEGE+HIGHSCHOOL+POOR+FAIR+GOOD+VGOOD+MNHPOOR+ANYLIMIT
      +HINCOME+MINCOME+LINCOME+NPOOR+INSURE,data=PersonLevel)
broom::tidy(glm_fit5)
broom::tidy(glm_fit5) %>% 
  dwplot(dot_args = list(size = 2, color = "black"),
         whisker_args = list(color = "black"),
         vline=geom_vline(xintercept=0,colour="grey50",linetype=2))

# Función para calcular heterocedasticidad corregida
SummaryML<-function(coefficients,n,k){ 
  # coeficientes=resultados dados por glm, 
  # n=número de observaciones, k=regresores incluyindo intercepto;
  Estimate=coefficients[,1]
  Std_Error=coefficients[,2]*sqrt((n-k)/n)
  t_Ratio=coefficients[,3]*sqrt(n/(n-k))
  coeffML=cbind(Estimate,Std_Error,t_Ratio)
  coeffML}
# Matriz de covarianza con heteroscedasticidad corregida; 
require(car)#install.packages("car"); ?hccm
cov_corrected=hccm(...,"hc1") #insertar Heckman - Etapa 2
var_corrected=diag(cov_corrected)
StdErr_corrected=sqrt(var_corrected)
tRatio_corrected=model2S$coefficients/StdErr_corrected
(summary_corrected=
  cbind(model2S$coefficients,StdErr_corrected,tRatio_corrected))
# Crear tabla con heterocedasticidad corregida ;
SummaryML(summary_corrected,2000,23)