1
\$\begingroup\$

I wrote a relatively simple R script to run a bunch of regressions. To be precise, it is a hurdle negative binomial regression (fancy stuff for count data). To print it out easily to Latex readable format, I store the results in a list.

I run five regressions in a row: one baseline regression, two specifications that add one variable at a time to the baseline, and finally two specifications that add one variable at a time to the four regressions before. The first two and the last two additional variables have the same beginning ("A" and "B" resp.) while the variables of the additional specifications have the same ending ("xyz" and "abc", resp.). To accomplish this, I construct the regression formula in a loop.

I also compare these regression results to the results of an ordinary negative binomial regression using Vuong's test (more fancy stuff on count data) right thereafter.

# Load packages
require(MASS) # for NegBin regression
require(pscl) # for Hurdle NegBin regression

# Load example data
Master <- read.csv("http://www.ats.ucla.edu/stat/data/fish.csv")
Master <- within(Master, {
  nofish <- factor(nofish)
  camper <- factor(camper)
  A_xyz <- persons
  A_abc <- persons*2
  B_xyz <- livebait
  B_abc <- livebait*2
})

# Define variables that are used throughout the script
status <- "without"
var_control <- c("child", "camper")
var_A <- c("A_xyz", "A_abc")
var_B <- c("B_xyz", "B_abc")

estimation_hurdle <- estimation_negbin <- list()

# Define a file to write out the results of Vuong's test
sink(file=paste0("filepath/", status, "/filename.txt"))

# Run baseline regression
writeLines("baseline")
estimation_hurdle$baseline <- hurdle(as.formula(formula_hurdle),
                                     data=Master, dist="negbin", link="probit"
)
estimation_negbin$baseline <- glm.nb(as.formula(formula_hurdle),
                                     data=Master
)
vuong(estimation_hurdle$baseline, estimation_negbin$baseline)

# Run regression for A variables
for (var in var_A) {
  writeLines(paste0("\n\n", var))
  estimation_hurdle[[var]] <- hurdle(as.formula(paste(formula_hurdle, var, sep=" + ")),
                                     data=Master, dist="negbin", link="probit"
  )
  estimation_negbin[[var]] <- glm.nb(as.formula(paste(formula_hurdle, var, sep=" + ")),
                                     data=Master
  )
  vuong(estimation_hurdle[[var]], estimation_negbin[[var]])
}
rm(var)

# Run regression for B and corresponding A variables
for (i in 1:length(var_B)) {
  writeLines(paste0("\n\n", var_B[i]))
  estimation_hurdle[[var_B[i]]] <- hurdle(as.formula(paste(formula_hurdle, var_A[i], var_B[i], sep=" + ")),
                                            data=Master, dist="negbin", link="probit"
  )
  estimation_negbin[[var_B[i]]] <- hurdle(as.formula(paste(formula_hurdle, var_A[i], var_B[i], sep=" + ")),
                                            data=Master
  )
  vuong(estimation_hurdle[[var_B[i]]], estimation_negbin[[var_B[i]]])
}
rm(i)
sink()

# Additional stuff to write out the results which are not much interesting here

Any idea to improve this code and make it easier/understandable?

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.