hetprobit.R
In hetprobit: Heteroscedastic Probit Regression Models

### R code from vignette source 'hetprobit.Rnw'

###################################################
### code chunk number 1: preliminaries
###################################################
options(width = 70, scipen = 1, digits = 4, prompt = "R> ", continue = "+  ")
library("hetprobit")
library("glmx")
library("lmtest")
library("memisc")


if(file.exists("VoterTurnout-models.rda")){
  load("VoterTurnout-models.rda")
} else {
# -----------------------------------------------+
#   Auxiliary Model Fitting File to speed up     |
#   (repeated) pdf-compilation of vignette       |
#------------------------------------------------+

data("VoterTurnout", package = "hetprobit")
library("hetprobit")

## homoscedastic model replicated from Nagler(1991)
mn <- glm(vote ~ age + I(age^2) + south + govelection + 
  (education + I(education^2)) * closing, 
  data = VoterTurnout, family = binomial(link = "probit"))

## heteroscedastic probit model
## with same regressors in mean and scale submodel
m1 <- hetprobit(vote ~ age + I(age^2) + south + govelection + 
  (education + I(education^2)) * closing |
  age + I(age^2) + south + govelection + 
  (education + I(education^2)) * closing, 
  data = VoterTurnout)

## reduced heteroscedastic model 
## excluding interaction effects
m2 <- update(m1, . ~ . - (education + I(education^2)):closing  |
  . - (education + I(education^2)):closing)

## reduced model for illustration purpose only
m <- hetprobit(vote ~ age + south + govelection + education + closing |
  age + south + govelection + education + closing, 
  data = VoterTurnout)

## benchmark model
## estimated with hetglm() from package glmx
m0 <- hetglm(vote ~ age + south + govelection + education + closing |
  age + south + govelection + education + closing, 
  data = VoterTurnout, method = "BFGS", hessian = TRUE)


save(mn, m, m0, m1, m2, file = "VoterTurnout-models.rda")
}

d <- mn$data


###################################################
### code chunk number 2: nagler1991-homoscedastic-probit (eval = FALSE)
###################################################
## data("VoterTurnout", package = "hetprobit")
## library("hetprobit")
## mn <- glm(vote ~ age + I(age^2) + south + govelection + 
##   (education + I(education^2)) * closing, 
##   data = VoterTurnout, family = binomial(link = "probit"))
## summary(mn)


###################################################
### code chunk number 3: hetprobit.Rnw:228-229
###################################################
summary(mn)


###################################################
### code chunk number 4: homoscedastic-probit (eval = FALSE)
###################################################
## mn1 <- hetprobit(vote ~ age + I(age^2) + south + govelection + 
##   (education + I(education^2)) * closing | 1, 
##   data = VoterTurnout)


###################################################
### code chunk number 5: heteroscedastic-probit-full (eval = FALSE)
###################################################
## m1 <- hetprobit(vote ~ age + I(age^2) + south + govelection + 
##   (education + I(education^2)) * closing |
##   age + I(age^2) + south + govelection + 
##   (education + I(education^2)) * closing, 
##   data = VoterTurnout)
## summary(m1)


###################################################
### code chunk number 6: hetprobit.Rnw:253-254
###################################################
summary(m1)


###################################################
### code chunk number 7: heteroscedastic-probit-reduced (eval = FALSE)
###################################################
## m2 <- update(m1, . ~ . - (education + I(education^2)):closing  |
##   . - (education + I(education^2)):closing)
## summary(m2)


###################################################
### code chunk number 8: hetprobit.Rnw:261-262
###################################################
summary(m2)


###################################################
### code chunk number 9: model-selection
###################################################
AIC(mn, m1, m2)
BIC(mn, m1, m2)


###################################################
### code chunk number 10: lr-test
###################################################
library("lmtest")
lrtest(mn, m1, m2)


###################################################
### code chunk number 11: glmx (eval = FALSE)
###################################################
## library("glmx")
## m0 <- hetglm(vote ~ age + south + govelection + education + 
##   I(education^2) + closing |
##   age + south + govelection + education + I(education^2) + closing, 
##   data = VoterTurnout, method = "BFGS", hessian = TRUE)


###################################################
### code chunk number 12: hetprobit.Rnw:302-303
###################################################
summary(m0)


###################################################
### code chunk number 13: hetprobit (eval = FALSE)
###################################################
## m <- hetprobit(vote ~ age + south + govelection + education + closing |
##   age + south + govelection + education + closing, 
##   data = VoterTurnout)


###################################################
### code chunk number 14: mtable (eval = FALSE)
###################################################
## library("memisc")
## toLatex(mtable(m, summary.stats = c("Log-likelihood", "AIC", "BIC", "N")))


###################################################
### code chunk number 15: mtable-latex
###################################################
toLatex(mtable(m, summary.stats = c("Log-likelihood", "AIC", "BIC", "N")))