MaxLikPMP: Estimates Bayesian Partial M-Probit via Maximum Likelihood.

Description Usage Arguments Details Value References See Also Examples

Description

MaxLikPMP obtains MLE estimates for a partial m-probit (pmp).

Usage

1
2
3
4
MaxLikPMP(formula, data, prpslid, gntid, R, I = NULL, V = NULL,
  qweights = NULL, pweights = NULL, qvote = NULL, pvote = NULL,
  betastart = NULL, verbose = TRUE, method = "BFGS", control = list(),
  whichloglik = NULL, betahat = NULL)

Arguments

formula

(required) a formula object of the form b ~ x1 + ... xK. All x and b must be without missing values. Intercept is included automatically. For each proposal, the b vector must be identical across actors.

data

(required) a data.frame object that contains all data used in the 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 actor in the data.

R

(required) integer. the number of actors that have to agree to pass a proposal excluding the actors that have veto power.

I

integer. The number of actors. Only required if max(gntid)==1,

V

integer. The number of actors that have have to agree to pass a proposal. The function assumes that the first V actors 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 actors with veto power.

betastart

a vector of starting value vectors for β. It is advised to supply starting values.

verbose

logical value. Use TRUE for progress report, FALSE for no output.

method

string value. Which of optim()'s optimization methods to use?

control

passed on to optim().

whichloglik

if set to "comb" the function will not make any attempt to find a more efficient representation of the likelihood function (see notes below)

q/pweights

integer matrix of dimension M x 1. The vote weights sorted by gntid.

q/pvote

scalar. The threshold for the vote weights.

behahat

a matrix of size K x H of known K coefficient values for the model for which the H likelihood values are computed

Details

Let x_{ij} be a vector of length K that collects all observed covariates for one of M actors and a proposal j. Let y_{ij} be the unobserved vote choice of actor i for proposal j. Let b_j be the observed vector of length J that collects the observed decisions made by the M actors according to a q-voting 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.

Different to BayesPMP(.) this function can not handle proposal-specific voting rules, partially observed voting records or varying intercepts.

Unless betahat is supplied, the function optimizes the implied likelihood using optim(). If betahat is supplied it computes the likelihood given using betahat as coefficients.

To make the computation as efficient as possible, it attempts to take advantage of two special cases of the partial m-probit likelihood.

a) If max(gntid)==1 and p/qweights==NULL the function assumes that all covariate values are the same across actors for each proposal. In this case the likelihood is equivalent to a Bernoulli likelihood with a special (voting-rule depended) link function that Marbach (2016) refers to as 'B-link'.

b) If p/qweights==NULL, the likelihood is the product of Poisson’s Binomial distribution functions which can be efficiently calculated using the discrete Fourier transform of the distribution's characteristic function (Hong 2013).

If none of the special cases applies the likelihood is computed by summing across the probabilities of all potential voting coalition implied by the voting rule. Even for moderately sized committees this computation is very time- and memory-intensive. To force the function to not attempt to detect a special case set whichloglik to "comb".

Value

maxlikpmp object with the following slots:

The standard R functions vcov, coef, logLik and summary can be used on any maxlikpmp object.

References

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

Hong, Yili. 2013. On Computing the Distribution Function for the Poisson Binomial Distribution. Computational Statistics & Data Analysis 59, 41–51.

See Also

optim

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
 ## Not run: 
 # Example 1: Using Simulated Data # 
 ###########################

 require(plyr)

 set.seed(10)
 J <- 250 	# proposal 
 I <- 10 	# actors
 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 <- MaxLikPMP(formula=y.agg ~ x1, prpslid="prpslid", gntid="gntid", data=data, I=I, R=R)

 summary(m1)

 # Plot log-likelihood function
 require(lattice)

 beta0 <- beta1 <- seq(-2,2,length=30)
 beta <- expand.grid(beta0,beta1)

 logliks <- MaxLikPMP(formula=y.agg ~ x1, prpslid="prpslid", gntid="gntid",
			     data=data, I=I, R=R, betahat=as.matrix(beta))

 logliks <- data.frame(loglik = logliks * (-1), beta0=beta[,1], beta1=beta[,2])

 wireframe(loglik ~ beta0 + beta1, data=logliks, 
	    shade=TRUE, main="Log-Likelihood", scales=list(arrows=FALSE) , zlim=c(-500,0) )


 # Example 2: Using Simulated Data # 
 ###########################

 No variation across actors for a proposal. Uses B-Link. 

 require(plyr)

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

 # Simualte recorded votes
 beta <- c(0,0.4)
 X <- data.frame(x0=1,x1=rep(runif(J,-2,2), I) ) # changed compared to Ex.1
 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 
 m2 <- MaxLikPMP(formula=y.agg ~ x1, prpslid="prpslid", gntid="gntid", data=data[data$gntid==1,], I=I, R=R)

 summary(m2)

 
## End(Not run)

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