BayesPMP: Estimates Bayesian Partial M-Probit via Gibbs Sampling.

In sumtxt/consilium: Estimating multivariate probit models with partial observability in R

Description

BayesPMP obtains a set of draws from the posterior distribution of a (mixed) partial m-probit.

Usage

BayesPMP(formula, data, prpslid, gntid, groups = NULL, R, V = NULL,
  qweights = NULL, pweights = NULL, qvote = NULL, pvote = NULL,
  pslct = NULL, qslct = NULL, b0 = NULL, invB0 = NULL, e0 = 0.1,
  f0 = 0.1, chib95 = FALSE, betastart = NULL, omegastart = NULL,
  burnin = 500, ngibbs = 2000, verbose = 500, thin = 1, chains = 2,
  seed = 42, rate = 200, step = 0.05, adapt = TRUE, monitor = FALSE)

Arguments

formula	(required) a formula object of the form y | b ~ x1 + ... xK. All x and b must be without missing values. For each proposal, the b vector must be identical across members. y | is optional supplying known vote choices for a proposal. Only used if all votes for a subset of proposal are none-missing.
data	(required) a data.frame object that contains all data used in formula and the variables specified in prpslid and gntid.
prpslid	(required) the name for the consecutive numbered ([1, J]) integer variable that identifies each proposal uniquely in the data.
gntid	(required) the name for the consecutive numbered ([1, M]) integer variable that identifies each voting member in the data.
groups	the name for the integer variable in data that identifies groups of proposals for which a varying intercept should be estimated.
R	(required) the number of members that have to agree to pass a proposal excluding the members that have veto power. Can be a vector of length n_distinct(prpslid) or a single number.
V	the number of veto members that have to agree to pass a proposal. Can can be a vector of length n_distinct(prpslid) or a single number. The function assumes that the first V members in the data are the once with veto power. That is, after sorting the data for each proposal by gntid the first V entries are used as they belong to members with with veto power.
b0	the prior mean of length K+1. If b0=NULL the function assumes b0=0.
invB0	the prior precision matrix of dimension (K+1) \times (K+1). If invB0=NULL the function assumes invB0=solve(diag(K+1)*1000) where K is the number of coefficients implied by formula.
chib95	calculate the marginal likelihood using the method of Chib (1995)?
betastart	the list of starting value vectors for β.
burnin	the integer number of samples used as burn-in.
ngibbs	the integer number of samples from the posterior distribution.
verbose	integer number: use 0 for no output; use any positive natural number to report progress at each verbose-th iteration.
thin	the n-th posterior draw to be recorded. Must be a number that yields a positive integer when ngibbs/thin.
chains	integer number for the number of chains.
seed	integer number for the seeding value.
rate	schedule for the adaption if adapt=TRUE.
step	ε-parameter if adapt=TRUE.
adapt	if TRUE (default) the Gibbs sampler uses the more efficient algorithm 2 in Marbach (2016) otherwise algorithm 1.
monitor	if adapt=TRUE prints adaptation steps for each draw.
q/pweights	integer matrix of dimension M x T where T is the number of distinct voting weights. The vote weights sorted by gntid.
q/pvote	scalar or integer vector of length J. The threshold for the vote weights.
q/pslct	integer vector of length J. Each entry refers to the applicable column in the q/pweights for a particular proposal j.
e0/f01	the parameters for the inverse-gamma prior density if a varying intercept is included (groups is not NULL)

Details

Let x_{ij} be a vector of length K that collects all observed covariates for one of M members and a proposal j. Let y_{ij} be the unobserved vote choice of member i for proposal j. Let b_j be the observed vector of length J that collects the observed decisions made by the M members according to a q-rule with threshold R. The model takes the following form:

Prob(b_j=0) = Prob( ∑_{i=1}^M y_{ij} < R)

Prob(y_{ij}=0) = φ(x_{ij}β)

where φ() is the standard univariate normal distribution and β is the parameter vector of interest. The prior density for β is a multivariate normal with user-specified prior mean vector and precision matrix. The function here simulates draws from the posterior density of β.

Notice, that BayesPMP can handle other voting rules than a q-rules with proposal-invariant threshold R. The threshold and the number of members with veto powers are allowed to vary across proposals. Weighted voting rules are also supported.

If a varying intercept is included, it is assumed to be drawn from a normal density with a variance that has an inverse-gamma prior density.

See the references for the derivation of the posterior and a description of the Gibbs sampler. The default (adapt=TRUE) uses algorithm 2 in Marbach (2016). The function runs in C++.

The marginal likelihood (chib95=TRUE) can not (yet) be calculated if a) a varying intercept is included (groups is NULL), b) the voting rule varies across decisions or c) with a partially observed voting record.

See coda documentation for help analyzing posterior samples. It is important to assess the convergence of the chains.

Value

coda object of ngibbs/thin draws from the posterior distribution of the coefficients including the intercept. The coda object has 3 non-standard attributes: type which is set to consilium, formula which is set to formula used in the function syntax and timecode recording the Sys.time() when the function finished. If chib95=TRUE the marginal likelihood estimate is also included.

#' @seealso summary.mcmc, plot.mcmc and other coda-functions.

References

Marbach, Moritz. 2016. 'Analyzing Decision Records from Committees.” Working Paper.

Chib, Siddhartha. 1995. Marginal Likelihood from the Gibbs Output. Journal of the American Statistical Association 90(432), 1313–1321.

Examples

 ## Not run: 
 # Example 1: q-rule # 
 ###########################

	require(plyr)

 set.seed(10)
 J <- 250 	# proposal 
 I <- 10 	# members
 R <- 6		# majority threshold

 # Simualte roll call voting record 
 beta <- c(0,0.4)
 X <- data.frame(x0=1,x1=runif(J*I,-2,2))
 y <- rbinom(J*I, 1, pnorm(as.matrix(X) %*% beta))

 # Bundle data with IDs
 data <- data.frame(gntid=sort(rep(seq(1,I), J)), 
 		prpslid=rep(seq(1,J), I), 
 		y, X)
 
 # Generate decision record 
 data <- ddply(data, "prpslid" ,function(x) { 
 		x$y.agg <- as.numeric(sum(x$y) >= R)
 		return(x)
 		})
 
 # Estimate partial m-probit 
 m1 <- BayesPMP(formula=y.agg ~ x1,R=R, prpslid="prpslid", gntid="gntid", data=data)

 # Generate selected individual voting record 
 data$y.slct <- ifelse(data$prpslid <= 10, data$y, NA)

 # Estimate partial m-probit with observed votes 
 m2 <- BayesPMP(formula=y.slct | y.agg ~ x1,R=R, prpslid="prpslid", gntid="gntid", data=data)
 
 summary(m1)
 summary(m2)

 # Example 2: Weighted q-rule 
 #############################

 require(plyr)
 
 set.seed(10)
 J <- 250 	# proposal 
 I <- 10 	# members
 R <- 6		# majority threshold
 voteweights <- matrix(0L, I, 1)
 voteweights[1:3,1] <- 10
 qvote <- 20
 qslct <- rep(1, J)

 # Simualte roll call voting record 
 beta <- c(0,0.4)
 X <- data.frame(x0=1,x1=runif(J*I,-2,2))
 y <- rbinom(J*I, 1, pnorm(as.matrix(X) %*% beta))

 # Bundle data with IDs
 data <- data.frame(gntid=sort(rep(seq(1,I), J)), 
 		prpslid=rep(seq(1,J), I), 
 		y, X)
 
 # Generate decision record 
 data <- ddply(data, "prpslid" ,function(x) { 
 		x$y.agg <- as.numeric( sum(x$y) >= R & (x$y %*% voteweights) >= qvote )
 		return(x)
 		})
 
 # Estimate partial m-probit 
 m2 <- BayesPMP(formula=y.agg ~ x1,R=R, prpslid="prpslid", 
		gntid="gntid", data=data, qweights=voteweights, qvote=qvote, qslct=qslct)
 summary(m2)
 
## End(Not run)

sumtxt/consilium documentation built on May 30, 2019, 8:38 p.m.