pscl: Political Science Computational Laboratory, Stanford University

CRAN: http://cran.r-project.org/web/packages/pscl/index.html
URL: http://pscl.stanford.edu/
Vignettes:
- Regression Models for Count Data in R

> library(pscl)

Loading required package: MASS
Loading required package: lattice
Classes and Methods for R developed in the

Political Science Computational Laboratory

Department of Political Science

Stanford University

Simon Jackman

hurdle and zeroinfl functions by Achim Zeileis

> data("bioChemists")

バージョン: 1.4.9

関数名	概略
`AustralianElectionPolling`	Political opinion polls in Australia, 2004-07
`AustralianElections`	elections to Australian House of Representatives, 1949-2007
`EfronMorris`	Batting Averages for 18 major league baseball players, 1970
`RockTheVote`	Voter turnout experiment, using Rock The Vote ads
`UKHouseOfCommons`	1992 United Kingdom electoral returns
`absentee`	Absentee and Machine Ballots in Pennsylvania State Senate Races
`admit`	Applications to a Political Science PhD Program
`betaHPD`	compute and optionally plot beta HDRs
`bioChemists`	article production by graduate students in biochemistry Ph.D. programs
`ca2006`	California Congressional Districts in 2006
`computeMargins`	add information about voting outcomes to a rollcall object
`constrain.items`	constrain item parameters in analysis of roll call data
`constrain.legis`	constrain legislators' ideal points in analysis of roll call data
`convertCodes`	convert entries in a rollcall matrix to binary form
`dropRollCall`	drop user-specified elements from a rollcall object
`dropUnanimous`	drop unanimous votes from rollcall objects and matrices
`extractRollCallObject`	return the roll call object used in fitting an ideal model
`hitmiss`	Table of Actual Outcomes against Predicted Outcomes for discrete data models
`hurdle`	Hurdle Models for Count Data Regression
`hurdle.control`	Control Parameters for Hurdle Count Data Regression
`hurdletest`	Testing for the Presence of a Zero Hurdle
`ideal`	analysis of educational testing data and roll call data with IRT models, via Markov chain Monte Carlo methods
`idealToMCMC`	convert an object of class ideal to a coda MCMC object
`igamma`	inverse-Gamma distribution
`iraqVote`	U.S. Senate vote on the use of force against Iraq, 2002.
`nj07`	rollcall object, National Journal key votes of 2007
`ntable`	nicely formatted tables
`odTest`	likelihood ratio test for over-dispersion in count data
`pR2`	compute various pseudo-R2 measures
`partycodes`	political parties appearing in the U.S. Congress
`plot.ideal`	plots an ideal object
`plot.predict.ideal`	plot methods for predictions from ideal objects
`plot.seatsVotes`	plot seats-votes curves
`politicalInformation`	Interviewer ratings of respondent levels of political information
`postProcess`	remap MCMC output via affine transformations
`predict.hurdle`	Methods for hurdle Objects
`predict.ideal`	predicted probabilities from an ideal object
`predict.zeroinfl`	Methods for zeroinfl Objects
`predprob`	compute predicted probabilities from fitted models
`predprob.glm`	Predicted Probabilties for GLM Fits
`predprob.ideal`	predicted probabilities from fitting ideal to rollcall data
`presidentialElections`	elections for U.S. President, 1932-2012, by state
`prussian`	Prussian army horse kick data
`readKH`	read roll call data in Poole-Rosenthal KH format
`rollcall`	create an object of class rollcall
`s109`	rollcall object, 109th U.S. Senate (2005-06).
`sc9497`	votes from the United States Supreme Court, from 1994-1997
`seatsVotes`	A class for creating seats-votes curves
`simpi`	Monte Carlo estimate of pi (3.14159265...)
`state.info`	information about the American states needed for U.S. Congress
`summary.ideal`	summary of an ideal object
`summary.rollcall`	summarize a rollcall object
`tracex`	trace plot of MCMC iterates, posterior density of legislators' ideal points
`unionDensity`	cross national rates of trade union density
`vectorRepresentation`	convert roll call matrix to series of vectors
`vote92`	Reports of voting in the 1992 U.S. Presidential election.
`vuong`	Vuong's non-nested hypothesis test
`zeroinfl`	Zero-inflated Count Data Regression
`zeroinfl.control`	Control Parameters for Zero-inflated Count Data Regression

bioChemists

> data("bioChemists")
> head(bioChemists)

  art   fem     mar kid5  phd ment
1   0   Men Married    0 2.52    7
2   0 Women  Single    0 2.05    6
3   0 Women  Single    0 3.75    6
4   0   Men Married    1 1.18    3
5   0 Women  Single    0 3.75   26
6   0 Women Married    2 3.59    2

> table(bioChemists$art)


  0   1   2   3   4   5   6   7   8   9  10  11  12  16  19 
275 246 178  84  67  27  17  12   1   2   1   1   2   1   1

hurdle

Arguments

formula
data, subset, na.action
weights
offest
dist
zero.dist
link
control
model, y, x
...

> # 次の２つのモデルの結果は等しい
> # logit-poisson
> (fm_hp1 <- bioChemists %>% hurdle(art ~ ., data = .))


Call:
hurdle(formula = art ~ ., data = .)

Count model coefficients (truncated poisson with log link):
(Intercept)     femWomen   marMarried         kid5          phd  
    0.67114     -0.22858      0.09648     -0.14219     -0.01273  
       ment  
    0.01875  

Zero hurdle model coefficients (binomial with logit link):
(Intercept)     femWomen   marMarried         kid5          phd  
    0.23680     -0.25115      0.32623     -0.28525      0.02222  
       ment  
    0.08012

> # geometric-poisson
> (fm_hp2 <- bioChemists %>% hurdle(art ~ ., data = ., zero = "geometric"))


Call:
hurdle(formula = art ~ ., data = ., zero.dist = "geometric")

Count model coefficients (truncated poisson with log link):
(Intercept)     femWomen   marMarried         kid5          phd  
    0.67114     -0.22858      0.09648     -0.14219     -0.01273  
       ment  
    0.01875  

Zero hurdle model coefficients (censored geometric with log link):
(Intercept)     femWomen   marMarried         kid5          phd  
    0.23680     -0.25115      0.32623     -0.28525      0.02222  
       ment  
    0.08012

> # logit-negative binominal
> (fm_hnb1 <- bioChemists %>% hurdle(art ~ ., data = ., dist = "negbin"))


Call:
hurdle(formula = art ~ ., data = ., dist = "negbin")

Count model coefficients (truncated negbin with log link):
(Intercept)     femWomen   marMarried         kid5          phd  
   0.355125    -0.244672     0.103417    -0.153260    -0.002933  
       ment  
   0.023738  
Theta = 1.8285 

Zero hurdle model coefficients (binomial with logit link):
(Intercept)     femWomen   marMarried         kid5          phd  
    0.23680     -0.25115      0.32623     -0.28525      0.02222  
       ment  
    0.08012

> # negbin-negbin
> (fm_hnb2 <- bioChemists %>% hurdle(art ~ ., data = ., dist = "negbin", zero = "negbin"))


Call:
hurdle(formula = art ~ ., data = ., dist = "negbin", zero.dist = "negbin")

Count model coefficients (truncated negbin with log link):
(Intercept)     femWomen   marMarried         kid5          phd  
   0.355125    -0.244672     0.103417    -0.153260    -0.002933  
       ment  
   0.023738  
Theta = 1.8285 

Zero hurdle model coefficients (censored negbin with log link):
(Intercept)     femWomen   marMarried         kid5          phd  
    8.06941     -2.36987      2.86178     -2.39751      0.05429  
       ment  
    0.84595  
Theta = 0.0746

zeroinfl

ゼロ過剰なカウントデータに対する最尤推定法による回帰モデル

Arguments

formula
data, subset, na.action
weights
offset
dist
link
control
model, y, x
...

> # 過剰を考慮しないモデル
> (fm_pois <- bioChemists %>% glm(art ~ ., data = ., family = poisson))


Call:  glm(formula = art ~ ., family = poisson, data = .)

Coefficients:
(Intercept)     femWomen   marMarried         kid5          phd  
    0.30462     -0.22459      0.15524     -0.18488      0.01282  
       ment  
    0.02554  

Degrees of Freedom: 914 Total (i.e. Null);  909 Residual
Null Deviance:        1817 
Residual Deviance: 1634     AIC: 3314

> (fm_zip <- bioChemists %>%  zeroinfl(art ~ . | 1, data = .))


Call:
zeroinfl(formula = art ~ . | 1, data = .)

Count model coefficients (poisson with log link):
(Intercept)     femWomen   marMarried         kid5          phd  
   0.553995    -0.231609     0.131971    -0.170474     0.002526  
       ment  
   0.021543  

Zero-inflation model coefficients (binomial with logit link):
(Intercept)  
     -1.681

> (fm_zinb <- bioChemists %>% zeroinfl(art ~ . | 1, data = ., dist = "negbin"))


Call:
zeroinfl(formula = art ~ . | 1, data = ., dist = "negbin")

Count model coefficients (negbin with log link):
(Intercept)     femWomen   marMarried         kid5          phd  
    0.25615     -0.21642      0.15049     -0.17642      0.01527  
       ment  
    0.02908  
Theta = 2.2644 

Zero-inflation model coefficients (binomial with logit link):
(Intercept)  
     -11.95