Commit 56aaff91 authored by Christopher Lawrence's avatar Christopher Lawrence Committed by Andreas Tille

Import Debian changes 1.03.10-1

r-cran-pscl (1.03.10-1) unstable; urgency=low

  * New upstream release.
parents c5a187ee 0e0ba8b2
Package: pscl
Version: 1.03.6.1
Date: 2010-08-25
Version: 1.03.10
Date: 2011-03-27
Title: Political Science Computational Laboratory, Stanford University
Author: Simon Jackman, with contributions from Alex Tahk, Achim
Zeileis, Christina Maimone and Jim Fearon
Maintainer: Simon Jackman <jackman@stanford.edu>
Depends: R (>= 2.11.0), MASS, stats, mvtnorm, coda, gam
Depends: R (>= 2.10.0), MASS, stats, mvtnorm, coda, gam, vcd
Suggests: MCMCpack, car, lmtest, sandwich, zoo
Enhances: stats, MASS
Imports: lattice
Description: Bayesian analysis of item-response theory (IRT) models,
roll call analysis; computing highest density regions; maximum
likelihood estimation of zero-inflated and hurdle models for
......@@ -18,6 +19,6 @@ LazyLoad: true
LazyData: true
License: GPL-2
URL: http://pscl.stanford.edu/
Packaged: 2011-03-28 23:03:05 UTC; jackman
Repository: CRAN
Date/Publication: 2011-03-21 14:58:50
Packaged: 2011-03-21 12:57:49 UTC; ripley
Date/Publication: 2011-03-31 15:18:56
1.03.10 * pythag deprecated in Rmath.h, use system hypot instead (3/13/2011)
* warnings about memory etc only come on with verbose=TRUE (req by Stephen Jessee)
1.03.9 * ideal: small change in partyLoyalty (thanks to Chris Hanretty)
* ideal: reformat output of ideal to be 3-d arrays
* ideal: change default prior precision for item parameters to .04 (used to be .01)
* added UKHouseOfCommons data; Example 6.9 in BASS
1.03.8 * added an optional "at" argument to predprob() methods for
count data so that the counts at which the probabilities
are evaluated can be specified
1.03.7 * small bug in constrain.item (reported by Paul Johnson)
* change normalization option in ideal to generate posterior means with mean 0, sd 1
* do normalization over all dimensions
* typos in documentation for pseudo-R2 (thanks to Henrik Pärn)
1.03.6 * made gam dependency explicit
* change linear.hypothesis to linearHypothesis
......@@ -11,7 +28,6 @@
* added nj07
* added vote92
1.03.3 * improved offset handling in hurdle()/zeroinfl(): offsets in zero
model are now allowed and can be different from count model.
See ?hurdle/?zeroinfl for details.
......
......@@ -572,7 +572,7 @@ model.matrix.hurdle <- function(object, model = c("count", "zero"), ...) {
}
predict.hurdle <- function(object, newdata, type = c("response", "prob", "count", "zero"),
na.action = na.pass, ...)
na.action = na.pass, at = NULL, ...)
{
type <- match.arg(type)
......@@ -629,7 +629,7 @@ predict.hurdle <- function(object, newdata, type = c("response", "prob", "count"
else if(!is.null(object$model)) y <- model.response(object$model)
else stop("predicted probabilities cannot be computed for fits with y = FALSE and model = FALSE")
yUnique <- min(y):max(y)
yUnique <- if(is.null(at)) 0:max(y) else at
nUnique <- length(yUnique)
rval <- matrix(NA, nrow = length(mu), ncol = nUnique)
dimnames(rval) <- list(rownames(X), yUnique)
......
......@@ -124,7 +124,7 @@ ideal <- function(object,
## check to see how much information will need to be stored
numrec <- (maxiter-burnin)/thin+1
if (interactive() &
if (interactive() & verbose &
((store.item)&&((n+m)*d*numrec>2000000))
||
((!store.item)&&((n*d*numrec)>2000000))
......@@ -138,7 +138,7 @@ ideal <- function(object,
stop("User terminated execution of ideal.")
}
if (interactive() & numrec>1000) {
if (interactive() & verbose & numrec>1000) {
ans <- readline(paste("You are attempting to save ",numrec," iterations. This\n",
"could result in a very large object and cause memory problems.\n",
"Do you want to continue with the current call to ideal? (y/n): ",
......@@ -222,8 +222,8 @@ ideal <- function(object,
else{
if(verbose)
cat("no prior precisions supplied for item parameters,\n",
"setting to default of .01\n")
bpv <- matrix(.01,m,d+1)
"setting to default of .04\n")
bpv <- matrix(.04,m,d+1)
}
if (((nrow(xp) != n)||(ncol(xp) != d)) || ((nrow(xpv)!=n)||(ncol(xpv)!=d))) {
......@@ -256,14 +256,14 @@ ideal <- function(object,
}
if(is.null(bpv)){
if(verbose)
cat("setting prior precisions for item parameters to all 0.01\n")
bpv <- matrix(0.01,m,d+1)
cat("setting prior precisions for item parameters to all 0.04\n")
bpv <- matrix(0.04,m,d+1)
}
xp <- as.vector(t(xp))
xpv <- as.vector(t(xpv))
bp <- as.vector(t(bp))
bpv <- as.vector(t(bpv))
xp <- as.vector(xp)
xpv <- as.vector(xpv)
bp <- as.vector(bp)
bpv <- as.vector(bpv)
################################################################
## check for start values - create if not supplied
......@@ -334,15 +334,25 @@ ideal <- function(object,
## report to user
if(verbose){
cat("using the following start values for ideal points (summary follows):\n")
print(summary(xstart))
cat("using the following start values for item parameters (summary follows):\n")
print(summary(bstart))
if(n<501){
cat("using the following start values for ideal points:\n")
print(xstart)
} else {
cat("using the following start values for ideal points (summary follows):\n")
print(summary(xstart))
}
if(m<501){
cat("using the following start values for item parameters:\n")
print(bstart)
} else {
cat("using the following start values for item parameters (summary follows):\n")
print(summary(bstart))
}
}
xstart <- as.vector(t(xstart))
bstart <- as.vector(t(bstart))
xstart <- as.vector(xstart)
bstart <- as.vector(bstart)
options(warn=0)
......@@ -351,7 +361,7 @@ ideal <- function(object,
##############################################################
yToC <- ifelse(is.na(v), 9, v)
yToC <- as.vector(t(yToC))
yToC <- as.vector(yToC)
cat("\nStarting MCMC Iterations...\n")
## ############################################
......@@ -420,58 +430,51 @@ ideal <- function(object,
cat("\n")
## parse returns from C job
## parse returns from C
xbar <- NULL
betabar <- NULL
if (!usefile) {
x <- output$xoutput
x <- matrix(x,nrow=numrec,byrow=T)
itervec <- seq(burnin,maxiter,by=thin)
x <- cbind(itervec,x)
rownames(x) <- x[,1]
colnames(x) <- gencolnames(legis.names,d)
keep <- itervec > burnin
## ideal points
print(output$xoutput[1:(n*d)])
x <- array(output$xoutput,
c(n,d,numrec))
## reshape to iteration first format
x <- aperm(x,c(3,1,2))
dimnames(x) <- list(itervec,
legis.names,
paste("D",1:d,sep=""))
if(verbose)
cat("...computing posterior means for ideal points...")
xbar <- getMean(keep,x)
if(verbose)
cat("done\n")
###############################################################
## item parameters
if (store.item) {
b <- output$boutput
b <- matrix(b,nrow=numrec,byrow=T)
b <- cbind(itervec,b)
rownames(b) <- b[,1]
colnames(b) <- gencolnames(vote.names,d,beta=T)
}
else {
b <- array(output$boutput,c(d+1,m,numrec)) ## parameters by votes by iters
dimnames(b) <- list(c(paste("Discrimination D",1:d,sep=""),
"Difficulty"),
vote.names,
itervec)
## reshape to iteration first format
b <- aperm(b,c(3,2,1)) ## iters by votes by parameters
if(verbose)
cat("...computing posterior means for item parameters...")
betabar <- getMean(keep,b)
if(verbose)
cat("done\n")
} else {
b <- NULL
}
}
else { ## output went to a file
} else { ## output went to a file
b <- x <- NULL
}
## compute some summary stats now, since we almost always use them
xbar <- betabar <- NULL
if(!is.null(x)){
if(verbose)
cat("MCMC sampling done, computing posterior means for ideal points...\n")
keep <- x[,1] > burnin
xbar <- apply(x[keep,-1],2,mean)
xbar <- matrix(xbar,n,d,byrow=TRUE)
mnames <- NULL
if(d>1){
for(j in 1:d)
mnames <- c(mnames,paste("Dimension",j))
}
dimnames(xbar) <- list(legis.names,mnames)
if(verbose)
cat("done\n")
}
if(store.item & !is.null(b)){
if(verbose)
cat("and for bill parameters...")
keep <- b[,1] > burnin
betabar <- apply(b[keep,-1],2,mean)
betabar <- matrix(betabar,m,d+1,byrow=TRUE)
if(verbose)
cat("done\n")
}
## wrap up for return to user
out <- list(n=n,m=m,d=d,
codes=codes,
......@@ -484,41 +487,17 @@ ideal <- function(object,
class(out) <- c("ideal")
## and, finally, if the user wanted meanzero
if(normalize)
if(normalize){
if(verbose)
cat("...normalizing output (post-processing)...")
out <- postProcess(out,
constraints="normalize")
if(verbose)
cat("done\n")
}
return(out)
}
gencolnames <- function(name, d, beta=FALSE) {
if(d>1){ ## more than one dimension?
dname <- NULL
for(i in 1:d){
dname <- c(dname,paste(name,"d",i,sep=""))
}
if(beta)
dname <- c(dname,paste(name,"Difficulty",sep=""))
dname <- matrix(dname,ncol=length(name),byrow=TRUE)
dname <- as.vector(dname)
dname <- c("Iteration",dname)
}
else {
if(beta){
if(beta)
dname <- c(name,paste(name,"Intercept",sep=""))
dname <- matrix(dname,ncol=length(name),byrow=T)
dname <- as.vector(dname)
dname <- c("Iteration",dname)
}
else {
dname <- c("Iteration",name)
}
}
dname
}
x.startvalues <- function(x,d,scale=TRUE,constraint=NULL,verbose=FALSE){
if(verbose)
cat("will use eigen-decomposition method to get start values for ideal points...")
......@@ -546,7 +525,7 @@ x.startvalues <- function(x,d,scale=TRUE,constraint=NULL,verbose=FALSE){
}
if(verbose)
cat("done\n")
v
return(v)
}
probit <- function(y,x){
......
......@@ -3,9 +3,10 @@
## check validity of a burnin number
## return logical of valid iters
checkBurnIn <- function(object, burnin) {
if (as.numeric(burnin)>=max(object$x[,1]))
stop("start must be less than the number of iterations")
return (object$x[,1] > burnin)
theIters <- as.numeric(dimnames(object$x)[[1]])
if (as.numeric(burnin)>max(theIters))
stop("burnin greater than number of iterations")
return (theIters > burnin)
}
checkD <- function(x,d) {
......@@ -18,28 +19,11 @@ checkCI <- function(conf.int) {
stop("conf.int must be between 0 and 1")
}
getDimX <- function(x,d,columns=TRUE) {
checkD(x,d)
px <- NULL
if(columns) {
px <- (x$x[,seq(from=d+1,to=ncol(x$x),by=x$d)])
colnames(px) <- x$legis.names
} else {
px <- (x$x[seq(from=d,to=nrow(x$x),by=x$d),])
rownames(px) <- x$legis.names
}
px
}
getDim <- function(x,d,dims,names,columns=TRUE) {
px <- NULL
if(columns) {
px <-(x[,seq(from=d,to=ncol(x),by=dims)])
colnames(px) <- names
}
else {
px <- (x[seq(from=d,to=nrow(x),by=dims),])
rownames(px) <- names
}
px
getMean <- function(keep,x){
xbar <- apply(x[keep,,,drop=FALSE],
c(2,3),
mean)
dimnames(xbar) <- list(dimnames(x)[[2]],
dimnames(x)[[3]])
return(xbar)
}
......@@ -43,11 +43,15 @@ plot1d <- function(x,
checkCI(conf.int) ## check that confidence interval is ok
q <- c((1-conf.int)/2, 1-((1-conf.int)/2)) ## quantiles from CI
xd <- getDimX(x,d) ## indicators for dimension
xm <- apply(xd[keep,],2,mean,na.rm=T) ## xbar
xm <- x$xbar ## xbar
indx <- order(xm) ## sort index
exispar <- par(no.readonly=T)
xq <- t(apply(xd[keep,],2,quantile,probs=q)) ## get CIs
exispar <- par(no.readonly=T)
myHPD <- function(x,prob){
tmp <- coda::as.mcmc(x)
return(coda::HPDinterval(tmp,prob))
}
xq <- t(apply(x$x[keep,,1],2,myHPD,prob=conf.int)) ## get HPDs
## names etc
cat(paste("Looking up legislator names and party affiliations\n"))
......@@ -60,13 +64,13 @@ plot1d <- function(x,
rm(tmpObject)
textLoc <- 1.05*min(xq) ## where to put x labels
if(showAllNames)
if(showAllNames){
par(mar=c(3,longName*.55,4,2)+0.1,
oma=rep(0,4))
else
} else {
par(mar=c(3,longName*.75,4,2)+0.1,
oma=rep(0,4))
}
## title string info
mainString <- paste("Ideal Points: ",
"Posterior Means and ",
......@@ -80,7 +84,8 @@ plot1d <- function(x,
x=1.02*range(xq),
xaxs="i",
xlab="",ylab="",
axes=F, type="n",
axes=FALSE,
type="n",
...)
mtext(mainString,side=3,line=3)
......@@ -107,7 +112,7 @@ plot1d <- function(x,
lines(y=c(i,i),x=xq[indx[i],],lwd=2)
if (is.null(party)){
points(y=i,x=xm[indx[i]],col="red",pch=19)
points(y=i,x=xm[indx[i]],col="red",pch=19,xpd=NULL)
}
else{
tbl <- table(party, exclude=NULL)
......@@ -116,7 +121,7 @@ plot1d <- function(x,
grp <- match(party[indx[i]],names(tbl))
points(y=i,
x=pt,
pch=19,col=cl[grp])
pch=19,col=cl[grp],xpd=NULL)
}
}
##par(ps=8)
......@@ -167,16 +172,13 @@ plot2d <- function(x,
if(d1==d2)
stop("can't do 2 dimensional summaries of the same dimension\n")
xd1 <- getDimX(x,d1)
xd2 <- getDimX(x,d2)
if(is.null(burnin)){ ## use x bar in ideal object
xm1 <- x$xbar[,d1]
xm2 <- x$xbar[,d2]
}
else{
xm1 <- apply(xd1[keep,],2,mean) ## posterior means
xm2 <- apply(xd2[keep,],2,mean)
xm1 <- apply(x$x[keep,,d1],2,mean) ## posterior means
xm2 <- apply(x$x[keep,,d2],2,mean)
}
if(oCP){
......@@ -186,10 +188,11 @@ plot2d <- function(x,
alphaBar <- x$betabar[,(x$d+1)]
}
else{
betaBar <- apply(x$beta[keep,-1],2,mean)
b1Bar <- betaBar[seq(from=d1,by=x$d+1,length=x$m)]
b2Bar <- betaBar[seq(from=d2,by=x$d+1,length=x$m)]
alphaBar <- betaBar[seq(from=x$d+1,by=x$d+1,length=x$m)]
bKeep <- x$beta[keep,,,drop=FALSE]
betaBar <- apply(bKeep,c(2,3),mean)
b1Bar <- betaBar[,d1]
b2Bar <- betaBar[,d2]
alphaBar <- betaBar[,x$d]
}
}
......@@ -204,6 +207,7 @@ plot2d <- function(x,
type="p",
xlab=paste("Dimension ",as.character(d1),sep=""),
ylab=paste("Dimension ",as.character(d2),sep=""),
xpd=NULL,
...)
if(oCP){
for(j in 1:x$m)
......@@ -212,8 +216,9 @@ plot2d <- function(x,
col=gray(.45))
}
}
else{
plot(x=xm1,y=xm2,
else{ ## we have party info
plot(x=xm1,
y=xm2,
main=mainString,
type="n",
xlab=paste("Dimension ",as.character(d1),sep=""),
......@@ -232,7 +237,8 @@ plot2d <- function(x,
thisParty <- party==names(tbl)[i]
points(y=xm2[thisParty],
x=xm1[thisParty],
pch=16,col=cl[i])
pch=16,col=cl[i],
xpd=NULL)
}
}
......@@ -315,7 +321,7 @@ tracex <- function(object,
keep <- checkBurnIn(object,eval(object$call$burnin,envir=.GlobalEnv))
else
keep <- checkBurnIn(object,burnin)
start <- object$x[keep,1][1]
start <- as.numeric(dimnames(object$x)[[1]])[keep][1]
## #######################################################
## one-dimensional stuff
......@@ -341,8 +347,8 @@ tracex <- function(object,
for (i in 1:nLegis){
meat <- object$x[keep,p[[i]]+d-1]
iter <- object$x[keep,1]
meat <- object$x[keep,p[[i]],1]
iter <- as.numeric(dimnames(object$x)[[1]])[keep]
par(mar=c(4, 4, 4, 2) + 0.1)
mainText <- plotName[i]
......@@ -411,8 +417,8 @@ tracex <- function(object,
col <- rainbow(nLegis) ## colors
meat <- list() ## container for iters to plot
for(i in 1:nLegis){
xTraces <- object$x[keep,p[[i]][1]]
yTraces <- object$x[keep,p[[i]][2]]
xTraces <- object$x[keep,p[[i]],d[1]]
yTraces <- object$x[keep,p[[i]],d[2]]
meat[[i]] <- list(x=xTraces,
y=yTraces,
col=col[i])
......@@ -434,8 +440,11 @@ tracex <- function(object,
axis(1,las=1)
axis(2,las=1)
lineFunc <- function(obj){
points(obj$x[1],obj$y[1],pch=16,col=obj$col)
lines(obj$x,obj$y,col=obj$col)
points(obj$x[1],obj$y[1],pch=1,col="black",cex=2)
npoints <- length(obj$x)
points(obj$x[npoints],obj$y[npoints],
pch=16,col="black",cex=2)
}
lapply(meat,lineFunc)
......
......@@ -4,26 +4,39 @@ postProcess <- function(object,
debug=FALSE){
if(class(object)!="ideal")
stop("postProcess only defined for objects of class ideal")
## process constraints
## not a list (i.e., usually "normalize" with d=1)
## not a list of normalizing constraint (i.e., usually "normalize=TRUE" with d=1)
if(!is.list(constraints)){
if(constraints=="normalize" & object$d>1)
stop("normalize option for postProcess only valid for one-dimesional models")
if(constraints=="normalize" & object$d==1){
## impose the mean zero, standard deviation one constraint
newObject <- normalizeIdeal(object)
if(constraints=="normalize"){
## get constraints needed for mean zero, standard deviation one restriction
tMat <- getNormalizingTransform(object) ## coefficients for linear map
if(debug){
cat("transformation matrix is:\n")
print(tMat)
}
newObject <- implementConstraints(object,tMat,debug)
}
}
if(is.list(constraints))
newObject <- postProcessAffine(object,constraints,debug)
return(newObject)
}
getNormalizingTransform <- function(object){
n <- object$n
m <- object$m
d <- object$d
offsets <- apply(object$xbar,2,mean)
s <- apply(object$xbar,2,sd)
coefs <- 1/s * diag(d) ## d by d
coefs <- rbind(coefs,-offsets/s) ## d+1 by d
return(coefs)
}
affineTrans <- function(x,target){
d <- dim(x)[2]
x0 <- cbind(x,1)
......@@ -39,67 +52,17 @@ affineTrans <- function(x,target){
}
foo <- solve(A)%*%b
foo <- matrix(foo,nrow=d+1)
foo <- cbind(foo,
c(rep(0,d),1))
foo
}
## normalize the ideal point estimates (d=1 only)
normalizeIdeal <- function(object){
d <- object$d
n <- object$n
m <- object$m
burnin <- eval(object$call$burnin)
theIters <- object$x[,1]
keep <- theIters > burnin
niter <- dim(object$x)[1]
haveBeta <- eval(object$call$store.item)
xMeans <- apply(object$x[,-1],1,mean)
xSd <- apply(object$x[,-1],1,sd)
xNew <- t(scale(t(object$x[,-1])))
if(haveBeta){
maxiter <- length(xMeans)
newBeta <- NA * object$beta
newBeta[,1] <- object$beta[,1]
d <- object$d
for(iter in 1:maxiter){
beta0 <- matrix(object$beta[iter,-1],
nrow=d+1,
byrow=FALSE)
beta0[2,] <- beta0[2,] - xMeans[iter]*beta0[1,] ## difficulty
beta0[1,] <- beta0[1,]*xSd[iter] ## discrimination
newBeta[iter,-1] <- as.vector(beta0)
}
}
newObject <- object
newObject$x <- cbind(object$x[,1],xNew)
dimnames(newObject$x) <- dimnames(object$x)
newObject$xbar <- matrix(apply(newObject$x[keep,-1],2,mean),
nrow=n,ncol=d,byrow=TRUE)
dimnames(newObject$xbar) <- dimnames(object$xbar)
if(haveBeta){
newObject$beta <- newBeta
newObject$betabar <- matrix(apply(newBeta[keep,-1],2,mean),
nrow=m,ncol=d+1,byrow=TRUE)
dimnames(newObject$betabar) <- dimnames(object$betabar)
}
newObject
return(foo)
}
postProcessAffine <- function(object,constraints,debug){
d <- object$d
n <- object$n
m <- object$m
burnin <- eval(object$call$burnin)
theIters <- object$x[,1]
keep <- theIters > burnin
niter <- dim(object$x)[1]
nSavedIters <- dim(object$x)[1]
theIters <- dimnames(object$x)[[1]]
keep <- checkBurnIn(object,
burnin=object$call$burnin)
nCon <- length(constraints)
if(nCon != (d+1)){
......@@ -115,7 +78,11 @@ postProcessAffine <- function(object,constraints,debug){
target <- matrix(NA,d+1,d)
for(i in 1:(d+1))
target[i,] <- constraints[[i]]
if(debug){
cat("target:\n")
print(target)
}
## form id vector, where are the named legislators in the ideal object?
legis.names <- dimnames(eval(object$call$object)$votes)[[1]]
if(is.null(legis.names)){
......@@ -142,65 +109,127 @@ postProcessAffine <- function(object,constraints,debug){
## initialize output objects
newX <- NA * object$x
newX[,1] <- theIters
dimnames(newX) <- dimnames(object$x)
haveBeta <- eval(object$call$store.item)
if(haveBeta){
cat("will also transform item/bill parameters\n")
newBeta <- NA * object$beta
newBeta[,1] <- object$beta[,1]
dimnames(newBeta) <- dimnames(object$beta)
}
## now loop over iterations
for(iter in 1:niter){
thisIter <- theIters[iter]
cat(paste("post-processing iteration",thisIter,"\n"))
x0 <- matrix(object$x[iter,-1],
ncol=d,
byrow=TRUE)