# racd10.R March 2012 for R version 2.15.0 # Create log file sink("racd10.Rout") rm(list=ls()) # ********** OVERVIEW OF racd10.R ********** # R Program # copyright C 2012 by A. Colin Cameron and Pravin K. Trivedi # used for "Regression Analyis of Count Data" SECOND EDITION # by A. Colin Cameron and Pravin K. Trivedi (2012) # Cambridge University Press # To run you need file # racd10data.dta # and R packages # foreign, sandwich # This R program reads in data and does Poisson regression for chapter 10 # (The Stata program racd10.do does more) # 10.5 EXAMPLE: DOCTOR VISITS # Poisson, NB2 and finite mixture of Poisson with 2 and 3 components # ********* DATA DESCRIPTION # The original data are from # P. Deb and P.K. Trivedi (2006a), "Specification and simulated # likelihood estimation of a nonnormal treatment-outcome model with # selection: application to health care utilization," # Econometrics Journal 9: 307-331. # See this article for more detailed discussion # ********* 10.5 DOCTOR VISITS: First column of Table 10.1 install.packages("foreign") library(foreign) data.ch10 <- read.dta(file = "racd10data.dta") data.ch10 <- subset(data.ch10, docvis<=70) # Selects <= 70 years data.ch10 <- subset(data.ch10, private*medicaid!=1) # Selects not both privte and Medicaid # Allows variables in database to be accessed simply by giving names attach(data.ch10) # Lists first six observations head(data.ch10) # Tabulate counts of doctor visits table(docvis) table(docvis) / nrow(data.ch10) # Variable list names(data.ch10) # Summary statistics summary(data.ch10) sapply(data.ch10,mean) sapply(data.ch10,sd) # Formula for the model estimated in this chapter - shortens the commands below formula.ch10model <- as.formula(docvis ~ private+medicaid+age+age2+educyr+actlim+totchr) # Poisson with default standard errors (variance is a multiple of the mean) model.poiss <- glm(formula.ch10model, family=poisson()) summary(model.poiss) ### First column of Table 10.1 # Poisson with robust standard errors (no assumption on variance) # install.packages("sandwich") library(sandwich) cov.robust <- vcovHC (model.poiss, type="HC0") se.robust <- sqrt(diag(cov.robust)) coeffs <- coef(model.poiss) t.robust <- coeffs / se.robust summary.poissrobust <- cbind(coeffs, se.robust, t.robust, pvalue=round(2*(1-pnorm(abs(coeffs/se.robust))),5), lower=coeffs-1.96*se.robust, upper=coeffs+1.96*se.robust) print(summary.poissrobust,digits=3) # close log file sink()