Title: | Discrete Choice (Binary, Poisson and Ordered) Models with Random Parameters |
---|---|
Description: | An implementation of simulated maximum likelihood method for the estimation of Binary (Probit and Logit), Ordered (Probit and Logit) and Poisson models with random parameters for cross-sectional and longitudinal data as presented in Sarrias (2016) <doi:10.18637/jss.v074.i10>. |
Authors: | Mauricio Sarrias [aut, cre] |
Maintainer: | Mauricio Sarrias <[email protected]> |
License: | GPL (>= 2) |
Version: | 0.3-6 |
Built: | 2025-02-27 03:16:28 UTC |
Source: | https://github.com/mauricio1986/rchoice |
Data from research by Long(1990) that analizes the scientist's level of publications.
data(Articles)
data(Articles)
A data frame with 915 observations on the following 6 variables:
art
Articles during last 3 years of Ph.D.,
fem
1 if female scientist; else 0,
mar
1 if married; else 0,
kid5
Number of children 5 or younger,
phd
Prestige of Ph.D. department,
ment
Articles by mentor during last 3 years,
Long, J. S. (1990). The origins of sex differences in science. Social Forces, 68(4), 1297-1316.
Long, J. S. (1997). Regression models for categorical and limited dependent variables (Vol. 7). Sage.
Long, J. S., & Freese, J. (2006). Regression models for categorical and limited dependent variables using Stata. Stata Press, College Station, TX.
data(Articles)
data(Articles)
In 1997 and 1989, the General Social Survey asked respondents to evaluate the following statement: "A working mother can establish just as warm and secure a relationship with her children as a mother who does not work".
data(Attitudes)
data(Attitudes)
A data frame with 2293 observations on the following 10 variables:
warm
1 = Strongly disagree, 2 = disagree, 3 = agree, 4 = strongly agree,
yr89
survey year: 1 = 1989; 0 = 1977,
male
1 = male; 0 = female,
white
1 = white; 0 = nonwhite,
age
age in years,
ed
years of education,
prst
occupational prestige,
Clogg, C. C., & Shihadeh, E. S. (1994). Statistical models for ordinal variables. Thousand Oaks, CA: Sage Publications.
Long, J. S. (1997). Regression models for categorical and limited dependent variables (Vol. 7). Sage.
Long, J. S., & Freese, J. (2006). Regression models for categorical and limited dependent variables using Stata. Stata Press, College Station, TX.
data(Attitudes)
data(Attitudes)
Computes the “bread” of the sandwich covariance matrix for a model of class Rchoice
## S3 method for class 'Rchoice' bread(x, ...)
## S3 method for class 'Rchoice' bread(x, ...)
x |
a fitted model of class |
... |
Other arguments when |
For more information see bread
from the package sandwich.
the covariance matrix times observations
Zeileis A (2006), Object-oriented Computation of Sandwich Estimators. Journal of Statistical Software, 16(9), 1–16.
## Probit model data("Workmroz") probit <- Rchoice(lfp ~ k5 + k618 + age + wc + hc + lwg + inc, data = Workmroz , family = binomial('probit')) summary(probit) library("sandwich") bread(probit)
## Probit model data("Workmroz") probit <- Rchoice(lfp ~ k5 + k618 + age + wc + hc + lwg + inc, data = Workmroz , family = binomial('probit')) summary(probit) library("sandwich") bread(probit)
Obtain the average marginal effects from hetprob
or ivpml
class models.
effect(object, ...)
effect(object, ...)
object |
an object of class |
... |
Additional arguments to be passed. |
Estimates of the average marginal effects computed as the average for each individual.
# Data library("AER") data("PSID1976") PSID1976$lfp <- as.numeric(PSID1976$participation == "yes") PSID1976$kids <- with(PSID1976, factor((youngkids + oldkids) > 0, levels = c(FALSE, TRUE), labels = c("no", "yes"))) PSID1976$finc <- PSID1976$fincome / 10000 # Average marginal effects for heteroskedastic Probit model labor_het <- hetprob(lfp ~ age + I(age^2) + finc + education + factor(kids) | factor(kids) + finc, data = PSID1976, link = "probit") eff_labor_het <- effect(labor_het) summary(eff_labor_het) # Average marginal effects for IV probit model # (nwincome is endogenous and heducation is the additional instrument) PSID1976$nwincome <- with(PSID1976, (fincome - hours * wage)/1000) fiml.probit <- ivpml(lfp ~ education + experience + I(experience^2) + age + youngkids + oldkids + nwincome | education + experience + I(experience^2) + age + youngkids + oldkids + heducation, data = PSID1976) summary(effect(fiml.probit)) summary(effect(fiml.probit, asf = FALSE))
# Data library("AER") data("PSID1976") PSID1976$lfp <- as.numeric(PSID1976$participation == "yes") PSID1976$kids <- with(PSID1976, factor((youngkids + oldkids) > 0, levels = c(FALSE, TRUE), labels = c("no", "yes"))) PSID1976$finc <- PSID1976$fincome / 10000 # Average marginal effects for heteroskedastic Probit model labor_het <- hetprob(lfp ~ age + I(age^2) + finc + education + factor(kids) | factor(kids) + finc, data = PSID1976, link = "probit") eff_labor_het <- effect(labor_het) summary(eff_labor_het) # Average marginal effects for IV probit model # (nwincome is endogenous and heducation is the additional instrument) PSID1976$nwincome <- with(PSID1976, (fincome - hours * wage)/1000) fiml.probit <- ivpml(lfp ~ education + experience + I(experience^2) + age + youngkids + oldkids + nwincome | education + experience + I(experience^2) + age + youngkids + oldkids + heducation, data = PSID1976) summary(effect(fiml.probit)) summary(effect(fiml.probit, asf = FALSE))
Obtain the average marginal effects from hetprob
class model.
## S3 method for class 'hetprob' effect(object, vcov = NULL, digits = max(3, getOption("digits") - 2), ...) ## S3 method for class 'effect.hetprob' summary(object, ...) ## S3 method for class 'effect.hetprob' print(x, ...) ## S3 method for class 'summary.effect.hetprob' print(x, digits = max(3, getOption("digits") - 3), ...)
## S3 method for class 'hetprob' effect(object, vcov = NULL, digits = max(3, getOption("digits") - 2), ...) ## S3 method for class 'effect.hetprob' summary(object, ...) ## S3 method for class 'effect.hetprob' print(x, ...) ## S3 method for class 'summary.effect.hetprob' print(x, digits = max(3, getOption("digits") - 3), ...)
object |
an object of class |
vcov |
an estimate of the asymptotic variance-covariance matrix of the parameters for a |
digits |
the number of digits. |
... |
further arguments.Ignored. |
x |
an object of class |
This function allows to obtain the average marginal effects (not the marginal effects at the mean). The standard errors are computed using Delta Method.
An object of class effect.heprob
.
Mauricio Sarrias.
# Data library("AER") data("PSID1976") PSID1976$lfp <- as.numeric(PSID1976$participation == "yes") PSID1976$kids <- with(PSID1976, factor((youngkids + oldkids) > 0, levels = c(FALSE, TRUE), labels = c("no", "yes"))) PSID1976$finc <- PSID1976$fincome / 10000 # Average marginal effects for heteroskedastic Probit model labor_het <- hetprob(lfp ~ age + I(age^2) + finc + education + factor(kids) | factor(kids) + finc, data = PSID1976, link = "probit") eff_labor_het <- effect(labor_het) summary(eff_labor_het)
# Data library("AER") data("PSID1976") PSID1976$lfp <- as.numeric(PSID1976$participation == "yes") PSID1976$kids <- with(PSID1976, factor((youngkids + oldkids) > 0, levels = c(FALSE, TRUE), labels = c("no", "yes"))) PSID1976$finc <- PSID1976$fincome / 10000 # Average marginal effects for heteroskedastic Probit model labor_het <- hetprob(lfp ~ age + I(age^2) + finc + education + factor(kids) | factor(kids) + finc, data = PSID1976, link = "probit") eff_labor_het <- effect(labor_het) summary(eff_labor_het)
Obtain the average marginal effects from ivpml
class model.
## S3 method for class 'ivpml' effect( object, vcov = NULL, asf = TRUE, digits = max(3, getOption("digits") - 2), ... ) ## S3 method for class 'effect.ivpml' summary(object, ...) ## S3 method for class 'effect.ivpml' print(x, ...) ## S3 method for class 'summary.effect.ivpml' print(x, digits = max(3, getOption("digits") - 3), ...)
## S3 method for class 'ivpml' effect( object, vcov = NULL, asf = TRUE, digits = max(3, getOption("digits") - 2), ... ) ## S3 method for class 'effect.ivpml' summary(object, ...) ## S3 method for class 'effect.ivpml' print(x, ...) ## S3 method for class 'summary.effect.ivpml' print(x, digits = max(3, getOption("digits") - 3), ...)
object |
an object of class |
vcov |
an estimate of the asymptotic variance-covariance matrix of the parameters for a |
asf |
if |
digits |
the number of digits. |
... |
further arguments.Ignored. |
x |
an object of class |
This function allows to obtain the average marginal effects (not the marginal effects at the mean). The standard errors are computed using Delta Method.
An object of class effect.ivpml
.
Mauricio Sarrias.
# Data library("AER") data("PSID1976") PSID1976$lfp <- as.numeric(PSID1976$participation == "yes") PSID1976$kids <- with(PSID1976, factor((youngkids + oldkids) > 0, levels = c(FALSE, TRUE), labels = c("no", "yes"))) # Average marginal effects for IV probit model # (nwincome is endogenous and heducation is the additional instrument) PSID1976$nwincome <- with(PSID1976, (fincome - hours * wage)/1000) fiml.probit <- ivpml(lfp ~ education + experience + I(experience^2) + age + youngkids + oldkids + nwincome | education + experience + I(experience^2) + age + youngkids + oldkids + heducation, data = PSID1976) summary(effect(fiml.probit)) summary(effect(fiml.probit, asf = FALSE))
# Data library("AER") data("PSID1976") PSID1976$lfp <- as.numeric(PSID1976$participation == "yes") PSID1976$kids <- with(PSID1976, factor((youngkids + oldkids) > 0, levels = c(FALSE, TRUE), labels = c("no", "yes"))) # Average marginal effects for IV probit model # (nwincome is endogenous and heducation is the additional instrument) PSID1976$nwincome <- with(PSID1976, (fincome - hours * wage)/1000) fiml.probit <- ivpml(lfp ~ education + experience + I(experience^2) + age + youngkids + oldkids + nwincome | education + experience + I(experience^2) + age + youngkids + oldkids + heducation, data = PSID1976) summary(effect(fiml.probit)) summary(effect(fiml.probit, asf = FALSE))
This a helper function to obtain the individuals' conditional estimate of the random parameters or compensating variations.
## S3 method for class 'Rchoice' effect(object, par = NULL, effect = c("cv", "ce"), wrt = NULL, ...)
## S3 method for class 'Rchoice' effect(object, par = NULL, effect = c("cv", "ce"), wrt = NULL, ...)
object |
an object of class |
par |
a string giving the name of the variable with random parameter, |
effect |
a string indicating what should be computed: the conditional expectation of the individual coefficients " |
wrt |
a string indicating respect to which variable the compensating variation should be computed, |
... |
further arguments. Ignored. |
A named list where “mean” contains the individuals' conditional mean for the random parameter or compensating variation, and where ‘sd.est’ contains their standard errors.
Greene, W. H. (2012). Econometric Analysis, Seventh Edition. Pearson Hall.
Train, K. (2009). Discrete Choice Methods with Simulation. Cambridge university press.
Rchoice
for the estimation of different discrete choice models with individual parameters.
# Poisson with random parameters data("Articles") poisson.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment, data = Articles, family = poisson, ranp = c(kid5 = "n", phd = "n", ment = "n"), R = 10) ## Get the individuals' conditional mean and their standard errors for ment bi.ment <- effect(poisson.ran, par = "ment", effect = "ce") summary(bi.ment$mean) summary(bi.ment$sd.est)
# Poisson with random parameters data("Articles") poisson.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment, data = Articles, family = poisson, ranp = c(kid5 = "n", phd = "n", ment = "n"), R = 10) ## Get the individuals' conditional mean and their standard errors for ment bi.ment <- effect(poisson.ran, par = "ment", effect = "ce") summary(bi.ment$mean) summary(bi.ment$sd.est)
It extracts the gradient for each observations evaluated at the estimated parameters for a model of class Rchoice
## S3 method for class 'Rchoice' estfun(x, ...)
## S3 method for class 'Rchoice' estfun(x, ...)
x |
a fitted model of class |
... |
Other arguments when |
For more information see estfun
from package sandwich.
the gradient matrix of dimension n times k
Zeileis A (2006), Object-oriented Computation of Sandwich Estimators. Journal of Statistical Software, 16(9), 1–16.
A generic function to collect coefficients and summary statistics from a effect.hetprob
object. It is used in mtable
## S3 method for class 'effect.hetprob' getSummary(obj, alpha = 0.05, ...)
## S3 method for class 'effect.hetprob' getSummary(obj, alpha = 0.05, ...)
obj |
an |
alpha |
level of the confidence intervals, |
... |
further arguments, |
For more details see package memisc.
A generic function to collect coefficients and summary statistics from a effect.ivpml
object. It is used in mtable
## S3 method for class 'effect.ivpml' getSummary(obj, alpha = 0.05, ...)
## S3 method for class 'effect.ivpml' getSummary(obj, alpha = 0.05, ...)
obj |
an |
alpha |
level of the confidence intervals, |
... |
further arguments, |
For more details see package memisc.
A generic function to collect coefficients and summary statistics from a hetprob
object. It is used in mtable
## S3 method for class 'hetprob' getSummary(obj, alpha = 0.05, ...)
## S3 method for class 'hetprob' getSummary(obj, alpha = 0.05, ...)
obj |
a |
alpha |
level of the confidence intervals, |
... |
further arguments, |
For more details see package memisc.
A generic function to collect coefficients and summary statistics from a ivpml
object. It is used in mtable
## S3 method for class 'ivpml' getSummary(obj, alpha = 0.05, ...)
## S3 method for class 'ivpml' getSummary(obj, alpha = 0.05, ...)
obj |
a |
alpha |
level of the confidence intervals, |
... |
further arguments, |
For more details see package memisc.
A generic function to collect coefficients and summary statistics from a Rchoice
object. It is used in mtable
## S3 method for class 'Rchoice' getSummary(obj, alpha = 0.05, ...)
## S3 method for class 'Rchoice' getSummary(obj, alpha = 0.05, ...)
obj |
a |
alpha |
level of the confidence intervals, |
... |
further arguments, |
For more details see package memisc.
German Health Care Data, unbalanced panel.
data(Health)
data(Health)
A data frame with 27326 observations on the following 27 variables:
id
person identification number
female
female =1, male =0
year
calendar year of the observation
age
age in years
hsat
health satisfaction, 0 (low),...,10 (high)
handdum
handicapped = 1, 0 otherwise
handper
degree of handicap in percent; 0,100
hhinc
household nominal monthly net income in German marks
hhkids
children under age 16 in the household = 1; otherwise = 0
educ
years of schooling
married
married =1, otherwise = 0
haupts
highest schooling degree is Hauptschul degree = 1; otherwise = 0
reals
highest schooling degree is Realschul degree = 1, otherwise = 0
fachhs
highest schooling degree is Polytechical degree = 1; otherwise = 0
abitur
highest schooling degree is Abitur = 1; otherwise = 0
univ
highest schooling degree is university degree =1; otherwise = 0
working
employed =1; otherwise = 0
bluec
blue-collar employee = 1; otherwise = 0
whitec
white-collar employeee =1; otherwise = 0
self
self-employed = 1; otherwise = 0
beamt
civil servant = 1; otherwise = 0
docvis
number of doctor visits in last three months
hospvis
number of hospital visits in last calendar year
public
insured in public health =1; otherwise = 0
addon
insured by add-on insurance =1; otherwise = 0
hsat2
40 observations on hsat recorded between 6 and 7 were changed to 7
newhsat
recording of hsat, (0-2) = 0, (3-5)=1, (6-8)=2, (9)=3 (10)=4
Riphahn, R. T., Wambach, A., & Million, A. (2003). Incentive effects in the demand for health care: a bivariate panel count data estimation. Journal of applied econometrics, 18(4), 387-405.
Greene, W. H. (2003). Econometric analysis. Pearson Education India.
data(Health)
data(Health)
Estimation of binary dependent variables, either probit or logit, with heteroskedastic error terms for cross-sectional dataset.
hetprob(formula, data, link = c("probit", "logit"), ...) ## S3 method for class 'hetprob' terms(x, ...) ## S3 method for class 'hetprob' model.matrix(object, ...) ## S3 method for class 'hetprob' estfun(x, ...) ## S3 method for class 'hetprob' bread(x, ...) ## S3 method for class 'hetprob' vcov(object, eigentol = 1e-12, ...) ## S3 method for class 'hetprob' df.residual(object, ...) ## S3 method for class 'hetprob' coef(object, ...) ## S3 method for class 'hetprob' logLik(object, ...) ## S3 method for class 'hetprob' print(x, ...) ## S3 method for class 'hetprob' summary(object, eigentol = 1e-12, ...) ## S3 method for class 'summary.hetprob' print(x, digits = max(3, getOption("digits") - 2), ...) ## S3 method for class 'hetprob' predict(object, newdata = NULL, type = c("xb", "pr", "sigma"), ...)
hetprob(formula, data, link = c("probit", "logit"), ...) ## S3 method for class 'hetprob' terms(x, ...) ## S3 method for class 'hetprob' model.matrix(object, ...) ## S3 method for class 'hetprob' estfun(x, ...) ## S3 method for class 'hetprob' bread(x, ...) ## S3 method for class 'hetprob' vcov(object, eigentol = 1e-12, ...) ## S3 method for class 'hetprob' df.residual(object, ...) ## S3 method for class 'hetprob' coef(object, ...) ## S3 method for class 'hetprob' logLik(object, ...) ## S3 method for class 'hetprob' print(x, ...) ## S3 method for class 'hetprob' summary(object, eigentol = 1e-12, ...) ## S3 method for class 'summary.hetprob' print(x, digits = max(3, getOption("digits") - 2), ...) ## S3 method for class 'hetprob' predict(object, newdata = NULL, type = c("xb", "pr", "sigma"), ...)
formula |
a symbolic description of the model of the form |
data |
the data of class |
link |
the assumption of the distribution of the error term. It could be either |
... |
arguments passed to |
x , object
|
an object of class |
eigentol |
the standard errors are only calculated if the ratio of the smallest and largest eigenvalue of the Hessian matrix is less than |
digits |
the number of digits. |
newdata |
optionally, a data frame in which to look for variables with which to predict. |
type |
the type of prediction required. The default, |
The heterokedastic binary model for cross-sectional data has the following structure:
with
where is the latent (unobserved) dependent variable for individual
;
is a
vector of independent variables determining the latent variable
(
x
variables in formula
);
and is the error term distributed either normally or logistically with
and heterokedastic variance
.
The variance for each individual is modeled parametrically assuming that it depends on a
vector observed variables
(
z
in formula
), whereas is the vector of parameters associated with each variable.
It is important to emphasize that
does not include a constant, otherwise the parameters are not identified.
The models are estimated using the maxLik
function from maxLik
package using both
analytic gradient and hessian (if Hess = TRUE
). In particular, the log-likelihood function is:
An object of class “hetprob
”, a list elements:
logLik0 |
logLik for the homokedastic model, |
f1 |
the formula, |
mf |
the model framed used, |
call |
the matched call. |
Mauricio Sarrias.
Greene, W. H. (2012). Econometric Analysis. 7 edition. Prentice Hall.
# Estimate a heteroskedastic probit and logit model data("Health") het.probit <- hetprob(working ~ factor(female) + factor(year) + educ + age + I(age^2) | factor(female) + age + I(age^2), data = Health, link = "probit") summary(het.probit) het.logit <- hetprob(working ~ factor(female) + factor(year) + educ + age + I(age^2) | factor(female) + age + I(age^2), data = Health, link = "logit") summary(het.logit)
# Estimate a heteroskedastic probit and logit model data("Health") het.probit <- hetprob(working ~ factor(female) + factor(year) + educ + age + I(age^2) | factor(female) + age + I(age^2), data = Health, link = "probit") summary(het.probit) het.logit <- hetprob(working ~ factor(female) + factor(year) + educ + age + I(age^2) | factor(female) + age + I(age^2), data = Health, link = "logit") summary(het.logit)
Estimation of Probit model with one endogenous and continuous variable by Maximum Likelihood.
ivpml(formula, data, messages = TRUE, ...) ## S3 method for class 'ivpml' terms(x, ...) ## S3 method for class 'ivpml' model.matrix(object, ...) ## S3 method for class 'ivpml' estfun(x, ...) ## S3 method for class 'ivpml' bread(x, ...) ## S3 method for class 'ivpml' vcov(object, ...) ## S3 method for class 'ivpml' df.residual(object, ...) ## S3 method for class 'ivpml' coef(object, ...) ## S3 method for class 'ivpml' logLik(object, ...) ## S3 method for class 'ivpml' print(x, ...) ## S3 method for class 'ivpml' summary(object, eigentol = 1e-12, ...) ## S3 method for class 'summary.ivpml' print(x, digits = max(3, getOption("digits") - 2), ...) ## S3 method for class 'ivpml' predict(object, newdata = NULL, type = c("xb", "pr", "stdp"), asf = TRUE, ...)
ivpml(formula, data, messages = TRUE, ...) ## S3 method for class 'ivpml' terms(x, ...) ## S3 method for class 'ivpml' model.matrix(object, ...) ## S3 method for class 'ivpml' estfun(x, ...) ## S3 method for class 'ivpml' bread(x, ...) ## S3 method for class 'ivpml' vcov(object, ...) ## S3 method for class 'ivpml' df.residual(object, ...) ## S3 method for class 'ivpml' coef(object, ...) ## S3 method for class 'ivpml' logLik(object, ...) ## S3 method for class 'ivpml' print(x, ...) ## S3 method for class 'ivpml' summary(object, eigentol = 1e-12, ...) ## S3 method for class 'summary.ivpml' print(x, digits = max(3, getOption("digits") - 2), ...) ## S3 method for class 'ivpml' predict(object, newdata = NULL, type = c("xb", "pr", "stdp"), asf = TRUE, ...)
formula |
a symbolic description of the model of the form |
data |
the data of class |
messages |
if |
... |
arguments passed to |
x , object
|
an object of class |
eigentol |
the standard errors are only calculated if the ratio of the smallest and largest eigenvalue of the Hessian matrix is less than |
digits |
the number of digits. |
newdata |
optionally, a data frame in which to look for variables with which to predict. |
type |
the type of prediction required. The default, |
asf |
if |
The IV probit for cross-sectional data has the following structure:
with
where is the latent (unobserved) dependent variable for individual
;
is the endogenous continuous variable;
is the vector of exogenous variables
which also includes the instruments for
;
and
are normal jointly distributed.
The model is estimated using the maxLik
function from maxLik
package using
analytic gradient.
Mauricio Sarrias.
Greene, W. H. (2012). Econometric Analysis. 7 edition. Prentice Hall.
# Data library("AER") data("PSID1976") PSID1976$lfp <- as.numeric(PSID1976$participation == "yes") PSID1976$kids <- with(PSID1976, factor((youngkids + oldkids) > 0, levels = c(FALSE, TRUE), labels = c("no", "yes"))) # IV probit model by MLE # (nwincome is endogenous and heducation is the additional instrument) PSID1976$nwincome <- with(PSID1976, (fincome - hours * wage)/1000) fiml.probit <- ivpml(lfp ~ education + experience + I(experience^2) + age + youngkids + oldkids + nwincome | education + experience + I(experience^2) + age + youngkids + oldkids + heducation, data = PSID1976) summary(fiml.probit)
# Data library("AER") data("PSID1976") PSID1976$lfp <- as.numeric(PSID1976$participation == "yes") PSID1976$kids <- with(PSID1976, factor((youngkids + oldkids) > 0, levels = c(FALSE, TRUE), labels = c("no", "yes"))) # IV probit model by MLE # (nwincome is endogenous and heducation is the additional instrument) PSID1976$nwincome <- with(PSID1976, (fincome - hours * wage)/1000) fiml.probit <- ivpml(lfp ~ education + experience + I(experience^2) + age + youngkids + oldkids + nwincome | education + experience + I(experience^2) + age + youngkids + oldkids + heducation, data = PSID1976) summary(fiml.probit)
Plot the distribution of the conditional expectation of the random parameters or compensating variations for objects of class Rchoice
.
## S3 method for class 'Rchoice' plot( x, par = NULL, effect = c("ce", "cv"), wrt = NULL, type = c("density", "histogram"), adjust = 1, main = NULL, col = "indianred1", breaks = 10, ylab = NULL, xlab = NULL, ind = FALSE, id = NULL, ... )
## S3 method for class 'Rchoice' plot( x, par = NULL, effect = c("ce", "cv"), wrt = NULL, type = c("density", "histogram"), adjust = 1, main = NULL, col = "indianred1", breaks = 10, ylab = NULL, xlab = NULL, ind = FALSE, id = NULL, ... )
x |
a object of class |
par |
a string giving the name of the variable with random parameter, |
effect |
a string indicating what should be plotted: the conditional expectation of the individual coefficients " |
wrt |
a string indicating repect to which variable should be computed the compensating variation, |
type |
a string indicating the type of distribution: it can be a |
adjust |
bandwidth for the kernel density, |
main |
an overall title for the plot, |
col |
color for the graph, |
breaks |
number of breaks for the histrogram if |
ylab |
a title for the y axis, |
xlab |
a title for the x axis, |
ind |
a boolean. If |
id |
only relevant if |
... |
further arguments. Ignored. |
Mauricio Sarrias
Greene, W. H. (2012). Econometric analysis, Seventh Edition. Pearson Hall.
Train, K. (2009). Discrete choice methods with simulation. Cambridge university press.
Rchoice
for the estimation of different discrete choice models with individual parameters.
# Poisson with random parameters data("Articles") poisson.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment, data = Articles, family = poisson, ranp = c(kid5 = "n", phd = "n", ment = "n"), R = 10) ## Plot the distribution of the conditional mean for ment plot(poisson.ran, par = "ment", type = "density") ## Plot the conditional mean for the first 20 individuals plot(poisson.ran, par = "ment", ind = TRUE, id = 1:20, col = "blue") ## Plot the compensating variation with respect to fem plot(poisson.ran, par = "ment", effect = "cv", wrt = "fem", type = "histogram")
# Poisson with random parameters data("Articles") poisson.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment, data = Articles, family = poisson, ranp = c(kid5 = "n", phd = "n", ment = "n"), R = 10) ## Plot the distribution of the conditional mean for ment plot(poisson.ran, par = "ment", type = "density") ## Plot the conditional mean for the first 20 individuals plot(poisson.ran, par = "ment", ind = TRUE, id = 1:20, col = "blue") ## Plot the compensating variation with respect to fem plot(poisson.ran, par = "ment", effect = "cv", wrt = "fem", type = "histogram")
Estimation of discrete choice models such as Binary (logit and probit), Poisson and Ordered (logit and probit) model with random coefficients for cross-sectional and panel data using simulated maximum likelihood.
Rchoice( formula, data, subset, weights, na.action, family, start = NULL, ranp = NULL, R = 40, haltons = NA, seed = NULL, correlation = FALSE, panel = FALSE, index = NULL, mvar = NULL, print.init = FALSE, init.ran = 0.1, gradient = TRUE, ... ) ## S3 method for class 'Rchoice' terms(x, ...) ## S3 method for class 'Rchoice' model.matrix(object, ...) ## S3 method for class 'Rchoice' coef(object, ...) ## S3 method for class 'Rchoice' fitted(object, ...) ## S3 method for class 'Rchoice' residuals(object, ...) ## S3 method for class 'Rchoice' df.residual(object, ...) ## S3 method for class 'Rchoice' update(object, new, ...) ## S3 method for class 'Rchoice' logLik(object, ...) ## S3 method for class 'Rchoice' print( x, digits = max(3, getOption("digits") - 3), width = getOption("width"), ... ) ## S3 method for class 'Rchoice' summary(object, ...) ## S3 method for class 'summary.Rchoice' print( x, digits = max(3, getOption("digits") - 3), width = getOption("width"), ... )
Rchoice( formula, data, subset, weights, na.action, family, start = NULL, ranp = NULL, R = 40, haltons = NA, seed = NULL, correlation = FALSE, panel = FALSE, index = NULL, mvar = NULL, print.init = FALSE, init.ran = 0.1, gradient = TRUE, ... ) ## S3 method for class 'Rchoice' terms(x, ...) ## S3 method for class 'Rchoice' model.matrix(object, ...) ## S3 method for class 'Rchoice' coef(object, ...) ## S3 method for class 'Rchoice' fitted(object, ...) ## S3 method for class 'Rchoice' residuals(object, ...) ## S3 method for class 'Rchoice' df.residual(object, ...) ## S3 method for class 'Rchoice' update(object, new, ...) ## S3 method for class 'Rchoice' logLik(object, ...) ## S3 method for class 'Rchoice' print( x, digits = max(3, getOption("digits") - 3), width = getOption("width"), ... ) ## S3 method for class 'Rchoice' summary(object, ...) ## S3 method for class 'summary.Rchoice' print( x, digits = max(3, getOption("digits") - 3), width = getOption("width"), ... )
formula |
a symbolic description of the model to be estimated. The |
data |
the data. It may be a |
subset |
an optional vector specifying a subset of observations, |
weights |
an optional vector of weigths, |
na.action |
a function wich indicated what should happen when the data
contains |
family |
the distribution to be used. It might be |
start |
a vector of starting values, |
ranp |
a named vector whose names are the random parameters and values the distribution:
" |
R |
the number of draws if |
haltons |
only relevant if |
seed |
the seed for the pseudo-random draws. This is only relevant if |
correlation |
only relevant if |
panel |
if |
index |
a string indicating the ‘id’ for individuals in the data. This argument is not required if data is a |
mvar |
only valid if |
print.init |
if |
init.ran |
initial values for standard deviation of random parameters. Default is 0.1, |
gradient |
if |
... |
further arguments passed to |
x , object
|
and object of class |
new |
an updated formula for the update method, |
digits |
number of digits, |
width |
width, |
The models are estimated using the maxLik
function from maxLik
package.
If ranp
is not NULL
, the random parameter model is estimated.
A random parameter model or random coefficient models permits regression parameter to
vary across individuals according to some distribution. A fully parametric
random parameter model specifies the latent variable conditional on regressors
and given parameters
to have conditional density
where
are iid with density
. The density is assumed a priori by the user by the argument
ranp
. If the parameters are assumed to be normally distributed , then the random parameter are constructed as:
where and
is the r-th draw from standard normal distribution for individual
.
Once the model is specified by the argument family
, the model is estimated using
Simulated Maximum Likelihood (SML). The probabilities, given by , are simulated using
R
pseudo-draws if halton=NULL
or R
halton draws if halton = NA
. The user can also specified the primes and the number of dropped elements for the halton draws. For example, if the model consists of two random parameters, the user can specify haltons = list("prime" = c(2, 3), "drop" = c(11, 11))
.
A random parameter hierarchical model can be estimated by including heterogeneity in the mean of the random parameters:
Rchoice manages the variables in the hierarchical model by the formula
object: all the hierarchical variables () are included after the
|
symbol. The argument mvar
indicate which variables enter in each random parameter. See examples below
An object of class “Rchoice
”, a list elements:
coefficients |
the named vector of coefficients, |
family |
type of model, |
link |
distribution of the errors, |
logLik |
a set of values of the maximum likelihood procedure, |
mf |
the model framed used, |
formula |
the formula (a Formula object), |
time |
|
freq |
frequency of dependent variable, |
draws |
type of draws used, |
R.model |
|
R |
number of draws used, |
bi |
an array of dimension |
Qir |
matrix of dimension |
ranp |
vector indicating the variables with random parameters and their distribution, |
probabilities |
the fitted probabilities for each individuals, |
residuals |
the residuals, |
call |
the matched call. |
Mauricio Sarrias [email protected]
Greene, W. H. (2012). Econometric Analysis. 7 edition. Prentice Hall.
Train, K. (2009). Discrete Choice Methods with Simulation. Cambridge university press.
## Probit model data("Workmroz") probit <- Rchoice(lfp ~ k5 + k618 + age + wc + hc + lwg + inc, data = Workmroz, family = binomial('probit')) summary(probit) ## Poisson model data("Articles") poisson <- Rchoice(art ~ fem + mar + kid5 + phd + ment, data = Articles, family = poisson) summary(poisson) ## Ordered probit model data("Health") oprobit <- Rchoice(newhsat ~ age + educ + hhinc + married + hhkids, data = Health, family = ordinal('probit'), subset = year == 1988) summary(oprobit) ## Poisson Model with Random Parameters poisson.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment, data = Articles, family = poisson, ranp = c(kid5 = "n", phd = "n", ment = "n")) summary(poisson.ran) ## Poisson Model with Correlated Random Parameters poissonc.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment, data = Articles, ranp = c(kid5 = "n", phd = "n", ment = "n"), family = poisson, correlation = TRUE, R = 20) summary(poissonc.ran) ## Hierarchical Poisson Model poissonH.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment | fem + phd, data = Articles, ranp = c(kid5 = "n", phd = "n", ment = "n"), mvar = list(phd = c("fem"), ment = c("fem", "phd")), family = poisson, R = 10) summary(poissonH.ran) ## Ordered Probit Model with Random Effects and Random Parameters Health$linc <- log(Health$hhinc) oprobit.ran <- Rchoice(newhsat ~ age + educ + married + hhkids + linc, data = Health[1:2000, ], family = ordinal('probit'), ranp = c(constant = "n", hhkids = "n", linc = "n"), panel = TRUE, index = "id", R = 10, print.init = TRUE) summary(oprobit.ran)
## Probit model data("Workmroz") probit <- Rchoice(lfp ~ k5 + k618 + age + wc + hc + lwg + inc, data = Workmroz, family = binomial('probit')) summary(probit) ## Poisson model data("Articles") poisson <- Rchoice(art ~ fem + mar + kid5 + phd + ment, data = Articles, family = poisson) summary(poisson) ## Ordered probit model data("Health") oprobit <- Rchoice(newhsat ~ age + educ + hhinc + married + hhkids, data = Health, family = ordinal('probit'), subset = year == 1988) summary(oprobit) ## Poisson Model with Random Parameters poisson.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment, data = Articles, family = poisson, ranp = c(kid5 = "n", phd = "n", ment = "n")) summary(poisson.ran) ## Poisson Model with Correlated Random Parameters poissonc.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment, data = Articles, ranp = c(kid5 = "n", phd = "n", ment = "n"), family = poisson, correlation = TRUE, R = 20) summary(poissonc.ran) ## Hierarchical Poisson Model poissonH.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment | fem + phd, data = Articles, ranp = c(kid5 = "n", phd = "n", ment = "n"), mvar = list(phd = c("fem"), ment = c("fem", "phd")), family = poisson, R = 10) summary(poissonH.ran) ## Ordered Probit Model with Random Effects and Random Parameters Health$linc <- log(Health$hhinc) oprobit.ran <- Rchoice(newhsat ~ age + educ + married + hhkids + linc, data = Health[1:2000, ], family = ordinal('probit'), ranp = c(constant = "n", hhkids = "n", linc = "n"), panel = TRUE, index = "id", R = 10, print.init = TRUE) summary(oprobit.ran)
Two kind of variables are used in models with individual heterogenetiy: the typical
variables that enter in the latent process and those variables that enter in the random
parameter (Hierarchical Model). rFormula
deal with this type of models using
suitable methods to extract the elements of the model.
rFormula(object) is.rFormula(object) ## S3 method for class 'rFormula' model.frame(formula, data, ..., lhs = NULL, rhs = NULL) ## S3 method for class 'rFormula' model.matrix(object, data, rhs = NULL, ...)
rFormula(object) is.rFormula(object) ## S3 method for class 'rFormula' model.frame(formula, data, ..., lhs = NULL, rhs = NULL) ## S3 method for class 'rFormula' model.matrix(object, data, rhs = NULL, ...)
object |
a formula form the |
formula |
a |
data |
a |
... |
further arguments. |
lhs |
see |
rhs |
see |
The vcov
method for Rchoice
objects extracts the covariance matrix of the coefficients or the random parameters. It also allows to get the standard errors for the variance-covariance matrix of the random parameters
## S3 method for class 'Rchoice' vcov( object, what = c("coefficient", "ranp"), type = c("cov", "cor", "sd"), se = FALSE, digits = max(3, getOption("digits") - 2), ... ) cov.Rchoice(x) cor.Rchoice(x) se.cov.Rchoice(x, sd = FALSE, digits = max(3, getOption("digits") - 2))
## S3 method for class 'Rchoice' vcov( object, what = c("coefficient", "ranp"), type = c("cov", "cor", "sd"), se = FALSE, digits = max(3, getOption("digits") - 2), ... ) cov.Rchoice(x) cor.Rchoice(x) se.cov.Rchoice(x, sd = FALSE, digits = max(3, getOption("digits") - 2))
object |
a fitted model of class |
what |
indicates which covariance matrix has to be extracted. The default is |
type |
if the model is estimated with random parameters, then this argument indicates what matrix should be returned. If |
se |
if |
digits |
number of digits, |
... |
further arguments |
x |
a fitted model of class |
sd |
if |
This new interface replaces the cor.Rchoice
, cov.Rchoice
and se.cov.Rchoice
functions which are deprecated.
Rchoice
for the estimation of discrete choice models with random parameters.
Data extracted by Mroz(1987) from the 1976 Panel Study of Income Dynacmis. The sample consists of 753 white, married women between the ages of 30 and 60.
data(Workmroz)
data(Workmroz)
A data frame with 753 observations on the following 9 variables:
lfp
1 if wife is in the paid labor force; else 0,
k5
Number of children ages 5 and younger,
k618
Number of children ages 6 to 18,
age
Wife's age in years,
wc
1 if wife attended college; else 0,
hc
1 if husband attended college; else 0,
lwg
Log of wife's estimated wage rate,
inc
Family income excluding wife's wage,
linc
Log of Family income excluding wife's wage,
Mroz, T. A. (1987). The sensitivity of an empirical model of married women's hours of work to economic and statistical assumptions. Econometrica, 55(4), 765-799
data(Workmroz)
data(Workmroz)