[r-cran-zelig] 12/102: Import Upstream version 2.4-1
Andreas Tille
tille at debian.org
Sun Jan 8 16:58:09 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-zelig.
commit 146fa76b968275574430f2ebf21a7daa447a2418
Author: Andreas Tille <tille at debian.org>
Date: Sun Jan 8 09:39:00 2017 +0100
Import Upstream version 2.4-1
---
DESCRIPTION | 10 +-
NAMESPACE | 10 +-
R/MCMC/class.ind.R | 11 ++
R/MCMC/qi.MCMCZelig.R | 242 ++++++++++++++++++++++--
R/MCMC/zelig2MCMC.R | 212 +++++++++++++++++----
R/MCMC/zelig3MCMC.R | 56 +++++-
R/class.ind.R | 11 ++
R/gee/param.gee.R | 34 ++++
R/gee/qi.gee.R | 69 +++++++
R/gee/vcov.gee.R | 3 +
R/gee/zelig2gamma.gee.R | 8 +
R/gee/zelig2logit.gee.R | 8 +
R/gee/zelig2normal.gee.R | 8 +
R/gee/zelig2poisson.gee.R | 8 +
R/load.first.R | 10 +-
R/mcmcei.R | 15 ++
R/param.MCMCZelig.R | 14 ++
R/param.relogit.R | 134 ++------------
R/print.relogit.R | 31 +---
R/print.relogit2.R | 7 +
R/print.summary.MCMCZelig.R | 13 ++
R/print.summary.glm.robust.R | 69 +------
R/print.summary.relogit.R | 14 ++
R/print.summary.relogit2.R | 6 +
R/qi.MCMCZelig.R | 409 +++++++++++++++++++++++++++++++++++++++++
R/qi.relogit.R | 32 ++--
R/qi.survreg.R | 109 ++++++-----
R/relogit.R | 123 +++++++++----
R/setx.relogit2.R | 7 +
R/sim.cond.R | 15 +-
R/sim.default.R | 6 +-
R/summarize.ei.R | 24 +++
R/summary.MCMCZelig.R | 19 ++
R/summary.relogit.R | 43 ++---
R/summary.relogit2.R | 21 +++
R/systemfit/callsystemfit.R | 8 +
R/systemfit/cmsystemfit.R | 19 ++
R/systemfit/param.multiple.R | 19 ++
R/systemfit/plot.zelig.2sls.R | 14 ++
R/systemfit/plot.zelig.3sls.R | 14 ++
R/systemfit/plot.zelig.sur.R | 14 ++
R/systemfit/plot.zelig.w2sls.R | 14 ++
R/systemfit/qi.multiple.R | 60 ++++++
R/systemfit/zelig22sls.R | 15 ++
R/systemfit/zelig23sls.R | 15 ++
R/systemfit/zelig2sur.R | 15 ++
R/systemfit/zelig2w2sls.R | 15 ++
R/vcov.relogit.R | 24 +--
R/zelig.R | 2 +-
R/zelig2MCMC.R | 348 +++++++++++++++++++++++++++++++++++
R/zelig2exp.R | 9 +-
R/zelig2logit.R | 4 +-
R/zelig2lognorm.R | 6 +-
R/zelig2negbin.R | 2 +-
R/zelig2probit.R | 2 +
R/zelig2relogit.R | 6 +-
R/zelig2weibull.R | 6 +-
R/zelig3MCMC.R | 136 ++++++++++++++
R/zelig3relogit.R | 53 ++++++
README | 22 +++
data/PErisk.txt | 63 +++++++
data/SupremeCourt.txt | 44 +++++
data/eidat.txt | 11 ++
data/newpainters.txt | 55 ++++++
data/swiss.txt | 48 +++++
data/tobin.txt | 21 +++
demo/00Index | 68 ++++---
demo/ei.dynamic.R | 55 ++++++
demo/ei.hier.R | 56 ++++++
demo/factor.bayes.R | 35 ++++
demo/factor.mix.R | 26 +++
demo/factor.ord.R | 28 +++
demo/irt1d.R | 32 ++++
demo/irtkd.R | 34 ++++
demo/logit.bayes.R | 57 ++++++
demo/match.R | 23 ++-
demo/mlogit.bayes.R | 58 ++++++
demo/normal.R | 24 +--
demo/normal.bayes.R | 58 ++++++
demo/oprobit.bayes.R | 67 +++++++
demo/poisson.bayes.R | 59 ++++++
demo/probit.bayes.R | 59 ++++++
demo/relogit.R | 39 +++-
demo/tobit.bayes.R | 59 ++++++
man/PErisk.Rd | 76 ++++++++
man/SupremeCourt.Rd | 31 ++++
man/eidat.Rd | 19 ++
man/newpainters.Rd | 40 ++++
man/swiss.Rd | 51 +++++
man/tobin.Rd | 32 ++++
90 files changed, 3489 insertions(+), 492 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 69260ec..1693876 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,13 +1,13 @@
Package: Zelig
-Version: 2.2-2
-Date: 2005-7-13
+Version: 2.4-1
+Date: 2005-08-15
Title: Zelig: Everyone's Statistical Software
Author: Kosuke Imai <kimai at Princeton.Edu>,
Gary King <king at harvard.edu>,
Olivia Lau <olau at fas.harvard.edu>
Maintainer: Olivia Lau <olau at fas.harvard.edu>
-Depends: R (>= 1.9.1), MASS, boot,
-Suggests: VGAM, survival, sandwich
+Depends: R (>= 2.0.0), MASS, boot,
+Suggests: VGAM, survival, sandwich, zoo, MCMCpack, coda
Description: Zelig is an easy-to-use program that can estimate, and
help interpret the results of, an enormous range of
statistical models. It literally is ``everyone's statistical
@@ -23,4 +23,4 @@ Description: Zelig is an easy-to-use program that can estimate, and
translates them into quantities of direct interest.
License: GPL version 2 or newer
URL: http://gking.harvard.edu/zelig
-Packaged: Wed Jul 13 16:09:18 2005; olau
+Packaged: Tue Aug 16 11:45:05 2005; olau
diff --git a/NAMESPACE b/NAMESPACE
index d70f148..f5b3847 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -56,11 +56,15 @@ S3method(print, BetaReg)
S3method(print, names.relogit)
S3method(print, names.zelig)
S3method(print, relogit)
+S3method(print, relogit2)
+S3method(print, summary.relogit)
+S3method(print, summary.relogit2)
S3method(print, summary.MI)
S3method(print, summary.strata)
S3method(print, summary.zelig.strata)
S3method(print, summary.lm.robust)
S3method(print, summary.glm.robust)
+S3method(print, summary.MCMCZelig)
S3method(print, summary.zelig)
S3method(print, zelig)
S3method(repl, zelig)
@@ -73,14 +77,18 @@ S3method(setx, default)
S3method(setx, noX)
S3method(setx, MI)
S3method(setx, strata)
+S3method(setx, relogit2)
S3method(summarize, array)
+S3method(summarize, ei)
S3method(summary, BetaReg)
S3method(summary, MI)
S3method(summary, strata)
S3method(summary, relogit)
+S3method(summary, relogit2)
S3method(summary, setx.cond)
S3method(summary, setx)
-S3method(summary, zelig.strata)
+S3method(summary, zelig.strata)
+S3method(summary, MCMCZelig)
S3method(summary, zelig)
S3method(summary, lm.robust)
S3method(summary, glm.robust)
diff --git a/R/MCMC/class.ind.R b/R/MCMC/class.ind.R
new file mode 100644
index 0000000..9c11983
--- /dev/null
+++ b/R/MCMC/class.ind.R
@@ -0,0 +1,11 @@
+#borrow from library(nnet)
+class.ind<-function (cl, levels.data=NULL)
+{
+ n <- length(cl)
+ cl <- as.factor(cl)
+ levels(cl)<-levels.data
+ x <- matrix(0, n, length(levels(cl)))
+ x[(1:n) + n * (unclass(cl) - 1)] <- 1
+ dimnames(x) <- list(names(cl), levels(cl))
+ x
+}
diff --git a/R/MCMC/qi.MCMCZelig.R b/R/MCMC/qi.MCMCZelig.R
index 5ef9608..924d43c 100644
--- a/R/MCMC/qi.MCMCZelig.R
+++ b/R/MCMC/qi.MCMCZelig.R
@@ -1,12 +1,15 @@
+
qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
model <- object$zelig
qi <- list()
check <- FALSE
- if (model %in% c("MCMClogit", "MCMCprobit"))
+
+ if (model %in% c("MCMClogit", "MCMCprobit", "MCMCoprobit", "MCMCmnl"))
check <- TRUE
- if (model %in% c("MCMClogit","MCMCprobit", "MCMCregress", "MCMCpoisson")) {
+ if (model %in% c("MCMClogit","MCMCprobit", "MCMCregress",
+ "MCMCpoisson","MCMCtobit")) {
if (model == "MCMClogit") {
coef <- object$coefficients
@@ -15,7 +18,8 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
dimnames(pr) <- dimnames(ev) <- dimnames(eta)
ev <- 1/(1+exp(-eta))
for (i in 1:ncol(ev))
- pr[,i] <- as.character(rbinom(length(ev[,i]), 1, ev[,i]))
+ pr[,i] <- as.character(rbinom(length(ev[,i]), 1, ev[,i]))
+
qi$ev <- ev
qi$pr <- pr
qi.name <- list(ev = "Expected Values: E(Y|X)", pr="Predicted
@@ -37,10 +41,32 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
else if (model =="MCMCregress") {
coef <- object$coefficients[,1:(ncol(object$coefficients)-1)]
eta <- coef %*% t(x)
- pr <- ev <- matrix(NA, nrow = nrow(eta), ncol = ncol(eta))
- dimnames(pr) <- dimnames(ev) <- dimnames(eta)
+ ev <- matrix(NA, nrow = nrow(eta), ncol = ncol(eta))
+ dimnames(ev) <- dimnames(eta)
ev <- eta
+ qi$ev <- ev
+ qi.name <- list(ev = "Expected Values: E(Y|X)")
+ }
+ else if (model =="MCMCtobit") {
+ coef <- object$coefficients[,1:(ncol(object$coefficients)-1)]
+ sig2 <- object$coefficients[,ncol(object$coefficients)]
+ sig <- sqrt(sig2)
+ eta <- coef %*% t(x)
+ ev <- cev <- matrix(NA, nrow = nrow(eta), ncol = ncol(eta))
+ dimnames(cev) <- dimnames(ev) <- dimnames(eta)
+
+ L2 <- (object$above-eta)/sig
+ L1 <- (object$below-eta)/sig
+ #cev <- eta + sig*(dnorm(L1)-dnorm(L2))/(pnorm(L2)-pnorm(L1))
+ if (object$below==-Inf) temp1<-0
+ else temp1 <- pnorm(L1)*object$below
+ if (object$above==Inf) temp2<-0
+ else temp2 <- (1-pnorm(L2))*object$above
+
+ ev <- temp1+eta*(pnorm(L2)-pnorm(L1))+sig*(dnorm(L1)-dnorm(L2))+temp2
+
qi$ev <- ev
+
qi.name <- list(ev = "Expected Values: E(Y|X)")
}
else if (model == "MCMCpoisson") {
@@ -58,7 +84,6 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
}
-
if (!is.null(x1)) {
eta1 <- coef %*% t(x1)
@@ -94,6 +119,24 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
}
+ else if (model == "MCMCtobit") {
+ L2 <- (object$above-eta1)/sig
+ L1 <- (object$below-eta1)/sig
+ #cev <- eta + sig*(dnorm(L1)-dnorm(L2))/(pnorm(L2)-pnorm(L1))
+
+ if (object$below==-Inf) temp1<-0
+ else temp1 <- pnorm(L1)*object$below
+ if (object$above==Inf) temp2<-0
+ else temp2 <- (1-pnorm(L2))*object$above
+
+ ev1 <- temp1+eta*(pnorm(L2)-pnorm(L1))+sig*(dnorm(L1)-dnorm(L2))+temp2
+
+ fd <-ev1-ev
+
+ qi$fd <- fd
+
+ qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
+ }
else if (model == "MCMCpoisson") {
ev1 <- exp(eta1)
fd <-ev1-ev
@@ -122,9 +165,177 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
}
list(qi=qi, qi.name=qi.name)
-}
-
- else if (model == "MCMChierEI") {
+ }
+ else if ((model =="MCMCoprobit") || (model == "MCMCmnl")) {
+ if (model == "MCMCoprobit") {
+ library(stats)
+ p <- dim(model.matrix(object, data=eval(object$data)))[2]
+ coef <- object$coefficients[,1:p]
+ eta <- coef %*% t(x)
+
+ level <- ncol(object$coefficients)-p+2
+ gamma<-matrix(NA, nrow(object$coefficients), level+1)
+
+ gamma[,1] <- rep(-Inf, nrow(gamma))
+ gamma[,2] <- rep(0, nrow(gamma))
+ gamma[,ncol(gamma)]<-rep(Inf, nrow(gamma))
+
+ if (ncol(gamma)>3)
+ gamma[,3:(ncol(gamma)-1)] <- object$coefficients[,(p+1):ncol(object$coefficients)]
+
+
+ ev <- array(NA, c(nrow(eta), level, ncol(eta)))
+ pr <- matrix(NA, nrow(eta), ncol(eta))
+# dimnames(pr)[1] <- dimnames(ev)[1] <- dimnames(eta)[1]
+# dimnames(pr)[2] <- dimnames(ev)[3] <- dimnames(eta)[2]
+
+ for (j in 1:level)
+ ev[,j,] <- pnorm(gamma[,j+1]-eta) - pnorm(gamma[,j]-eta)
+
+ for (j in 1:nrow(pr)) {
+ mu <- eta[j,]
+ pr[j,]<-as.character(cut(mu, gamma[j,], labels=as.factor(1:level)))
+ }
+
+
+ # t(t(1:nrow(pr))%*%pr)
+ qi$ev <- ev
+ qi$pr <- pr
+ qi.name <- list(ev = "Expected Values: E(Y|X)", pr="Predicted
+ Values: Y|X")
+ }
+ else if (model == "MCMCmnl") {
+ library(stats)
+ resp <- model.response(model.frame(object))
+ level <- length(table(resp))
+
+ p <- dim(model.matrix(eval(object),data=eval(object$data)))[2]
+
+ coef <- object$coefficients
+
+
+ eta <- array(NA, c(nrow(coef),level, nrow(x)))
+
+ eta[,1,]<-matrix(0, dim(eta)[1],dim(eta)[3])
+ for (j in 2:level) {
+ ind <- (1:p)*(level-1)-(level-j)
+ eta[,j,]<- coef[,ind]%*%t(x)
+ }
+
+ eta<-exp(eta)
+
+ ev <- array(NA, c(nrow(coef), level, nrow(x)))
+ pr <- matrix(NA, nrow(coef), nrow(x))
+ #dimnames(pr)[1] <- dimnames(ev)[1] <- dimnames(eta)[1]
+
+ for (k in 1:nrow(x))
+ for (j in 1:level)
+ {
+ ev[,j,k] <- eta[,j,k]/rowSums(eta[,,k])
+ }
+
+ for (k in 1:nrow(x)) {
+ probs <- as.matrix(ev[,,k])
+ temp <- apply(probs, 1, FUN=rmultinom, n=1, size=1)
+ temp <- as.matrix(t(temp)%*%(1:nrow(temp)))
+ pr <- apply(temp,2,as.character)
+ }
+ qi$ev <- ev
+ qi$pr <- pr
+ qi.name <- list(ev = "Expected Values: E(Y|X)", pr="Predicted
+ Values: Y|X")
+ }
+
+
+ if (!is.null(x1)) {
+
+ if (model == "MCMCoprobit") {
+
+ eta1 <- coef %*% t(x1)
+
+ ev1 <- array(NA, c(nrow(eta), level, ncol(eta)))
+
+ for (j in 1:level)
+ ev1[,j,] <- pnorm(gamma[,j+1]-eta1) - pnorm(gamma[,j]-eta1)
+
+ rr <-ev1/ev
+ fd <-ev1-ev
+
+ qi$fd <- fd
+ qi$rr <- rr
+
+ qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
+
+ qi.name$rr <- "Risk Ratios: P(Y=1|X1)/P(Y=1|X)"
+
+ }
+ else if (model == "MCMCmnl") {
+
+ eta1 <- array(NA, c(nrow(coef),level, nrow(x1)))
+
+ eta1[,1,]<-matrix(0, dim(eta1)[1],dim(eta1)[3])
+ for (j in 2:level) {
+ ind <- (1:p)*(level-1)-(level-j)
+ eta1[,j,]<- coef[,ind]%*%t(x1)
+ }
+
+ eta1<-exp(eta1)
+
+ ev1 <- array(NA, c(nrow(eta1), level, nrow(x1)))
+
+ for (k in 1:nrow(x))
+ for (j in 1:level)
+ {
+ ev1[,j,k] <- eta1[,j,k]/rowSums(eta1[,,k])
+ }
+
+ rr <-ev1/ev
+ fd <-ev1-ev
+
+ qi$fd <- fd
+ qi$rr <- rr
+
+ qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
+
+ qi.name$rr <- "Risk Ratios: P(Y=1|X1)/P(Y=1|X)"
+
+ }
+ }
+
+ if (!is.null(y)) {
+ yvar <- matrix(rep(y, nrow(simpar)), nrow = nrow(simpar), byrow =
+ TRUE)
+
+ levels.names<-levels(as.factor(y))
+ yvar1<- pr1 <-array(NA, c(nrow(yvar), level, ncol(yvar)))
+
+
+ for (j in 1:nrow(yvar))
+ {
+ yvar1[j,,]<-t(class.ind(yvar[j,], levels.names))
+ if (check)
+ pr1[j,,]<-t(class.ind(as.integer(qi$pr[j,]), levels.names))
+ else
+ pr1[j,,]<-t(class.ind(qi$pr[j,],levels.names))
+ }
+
+ tmp.ev <- yvar1 - qi$ev
+ tmp.pr <- yvar1 - pr1
+ qi$ate.ev <- matrix(apply(tmp.ev, 2, rowMeans), nrow = nrow(simpar))
+ qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
+ if (model %in% c("MCMCoprobit", "MCMCmnl"))
+ {
+ qi$ate.pr <- matrix(apply(tmp.pr, 2, rowMeans), nrow = nrow(simpar))
+ qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
+ }
+ }
+
+ list(qi=qi, qi.name=qi.name)
+
+
+
+ }
+ else if (model == "MCMChierEI" || model == "MCMCdynamicEI") {
if (!any(class(x)=="cond")) stop("set 'cond=TRUE' in setx.\n")
else
{
@@ -164,8 +375,17 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
list(qi=qi, qi.name=qi.name)
}
-}
-
+ else if ( model %in% c("MCMCfactanal", "MCMCordfactanal",
+ "MCMCmixfactanal", "MCMCirt1d", "MCMCirtKd"))
+ {
+ stop("sim procedure not applicable since no explanatory
+ variables are involved.\n")
+ list(qi=qi)
+ }
+
+# return(list(qi=qi, qi.name=qi.name))
+
+}
diff --git a/R/MCMC/zelig2MCMC.R b/R/MCMC/zelig2MCMC.R
index 91ef2a4..134cab3 100644
--- a/R/MCMC/zelig2MCMC.R
+++ b/R/MCMC/zelig2MCMC.R
@@ -13,7 +13,7 @@ zelig2MCMChierEI <- function(formula, model, data, M, ...) {
}
mf$model <- mf$M <- NULL
- temp <- mcmcei(formula=formula, data=data, covar=covar)
+ temp <- mcmcei(formula=formula, data=data)
if ((any(temp<0)) || ((any(temp<1) && !any(temp==0) ) && any(temp>1)))
stop("data need to be either counts or proportions.\n")
@@ -43,6 +43,53 @@ zelig2MCMChierEI <- function(formula, model, data, M, ...) {
as.call(mf)
}
+zelig2MCMCdynamicEI <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 50000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 5000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf$model <- mf$M <- NULL
+ temp <- mcmcei(formula=formula, data=data)
+
+ if ((any(temp<0)) || ((any(temp<1) && !any(temp==0) ) && any(temp>1)))
+ stop("data need to be either counts or proportions.\n")
+ if (is.null(mf$N))
+ {
+ if (all(temp>=0)) #N is not needed
+ {
+ mf$r0 <- temp$r0
+ mf$r1 <- temp$r1
+ mf$c0 <- temp$c0
+ mf$c1 <- temp$c1
+ }
+ else stop("Needs total counts for inputs as porportion.\n")
+ }
+ else if (((length(mf$N)!= nrow(data)) && (length(mf$N)!=1)) || (any(mf$N<1)))
+ stop("N needs to have same length as the observations and be postive numbers\n.")
+ else if ((all(temp<1)) && (all(mf$N>1)))
+ {
+ mf$r0 <- round(temp$r0*mf$N)
+ mf$r1 <- mf$N-mf$r0
+ mf$c0 <- round(temp$c0*mf$N)
+ mf$c1 <- mf$N-mf$c0
+
+ }
+
+ mf[[1]] <- MCMCpack::MCMCdynamicEI
+ as.call(mf)
+}
+
+
+
zelig2MCMClogit <- function(formula, model, data, M, ...) {
require(MCMCpack)
@@ -119,68 +166,86 @@ zelig2MCMCpoisson <- function(formula, model, data, M, ...) {
as.call(mf)
}
-zelig2MCMCoprobit <- function(formula, model, data, M, ...) {
+zelig2MCMCtobit <- function(formula, model, data, M, ...) {
require(MCMCpack)
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- NULL
if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
- else mf$verbose <- 500
- mf[[1]] <- MCMCpack::MCMCoprobit
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf$model <- mf$M <- NULL
+ mf[[1]] <- MCMCpack::MCMCtobit
as.call(mf)
}
-
-
-
zelig2MCMCmnl <- function(formula, model, data, M, ...) {
require(MCMCpack)
+ require(stats)
mf <- match.call(expand.dots = TRUE)
mf$model <- mf$M <- NULL
if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
- else mf$verbose <- 500
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
mf[[1]] <- MCMCpack::MCMCmnl
as.call(mf)
}
-
-
-
-zelig2MCMCtobit <- function(formula, model, data, M, ...) {
+zelig2MCMCoprobit <- function(formula, model, data, M, ...) {
require(MCMCpack)
+ require(stats)
mf <- match.call(expand.dots = TRUE)
- if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
- else mf$verbose <- 500
mf$model <- mf$M <- NULL
- mf[[1]] <- MCMCpack::MCMCtobit
- as.call(mf)
-}
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
-zelig2MCMCdynamicEI <- function(formula, model, data, M, ...) {
- require(MCMCpack)
- mf <- match.call(expand.dots = TRUE)
- if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
- else mf$verbose <- 500
- mf$model <- mf$M <- NULL
- temp <- mcmcei(formula=formula, data=data, covar=covar)
- mf$r0 <- temp$r0
- mf$r1 <- temp$r1
- mf$c0 <- temp$c0
- mf$c1 <- temp$c1
- # mf$formula <- mf$data <- NULL
- mf[[1]] <- MCMCpack::MCMCdynamicEI
+ mf[[1]] <- MCMCpack::MCMCoprobit
as.call(mf)
}
+
+
+
zelig2MCMCfactanal <- function(formula, model, data, M, ...) {
require(MCMCpack)
mf <- match.call(expand.dots = TRUE)
- if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
- else mf$verbose <- 500
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 20000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ if (is.null(mf$factors)) mf$factors<-2
+ else if (mf$factors<2) stop("number of factors needs to be at
+ least 2")
mf$model <- mf$M <- NULL
- mf$x <- as.matrix(model.frame(formula, data=data, na.action=NULL))
+ mf$x <- as.matrix(model.response(model.frame(formula, data=data, na.action=NULL)))
mf[[1]] <- MCMCpack::MCMCfactanal
as.call(mf)
}
@@ -188,10 +253,21 @@ zelig2MCMCfactanal <- function(formula, model, data, M, ...) {
zelig2MCMCordfactanal <- function(formula, model, data, M, ...) {
require(MCMCpack)
mf <- match.call(expand.dots = TRUE)
- if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
- else mf$verbose <- 500
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 20000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+ if (is.null(mf$factors)) mf$factors<-2
+ else if (mf$factors<1) stop("number of factors needs to be at
+ least 1")
mf$model <- mf$M <- NULL
- mf$x <- as.matrix(model.frame(formula, data=data, na.action=NULL))
+ mf$x <- as.matrix(model.response(model.frame(formula, data=data, na.action=NULL)))
+
mf[[1]] <- MCMCpack::MCMCordfactanal
as.call(mf)
}
@@ -200,13 +276,73 @@ zelig2MCMCmixfactanal <- function(formula, model, data, M, ...) {
require(MCMCpack)
mf <- match.call(expand.dots = TRUE)
if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
- else mf$verbose <- 500
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ if (is.null(mf$factors)) mf$factors<-2
+ else if (mf$factors<1) stop("number of factors needs to be at
+ least 1")
mf$model <- mf$M <- NULL
+
var <- model.response(model.frame(formula, data=data,
na.action=NULL))
varnames <- colnames(var)
mf$x <- as.formula(paste("~", paste(varnames, collapse="+")))
mf[[1]] <- MCMCpack::MCMCmixfactanal
+
+ as.call(mf)
+}
+
+
+zelig2MCMCirt1d <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 20000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf$model <- mf$M <- NULL
+ mf$datamatrix <- as.matrix(model.response(model.frame(formula, data=data,
+ na.action=NULL)))
+ mf$datamatrix <- t(mf$datamatrix)
+ mf[[1]] <- MCMCpack::MCMCirt1d
+
+ as.call(mf)
+}
+
+
+zelig2MCMCirtKd <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+ if (is.null(mf$dimensions)) mf$dimensions <- 1
+
+ mf$model <- mf$M <- NULL
+ mf$datamatrix <- as.matrix(model.response(model.frame(formula, data=data,
+ na.action=NULL)))
+ mf$datamatrix <- t(mf$datamatrix)
+ mf[[1]] <- MCMCpack::MCMCirtKd
+
as.call(mf)
}
diff --git a/R/MCMC/zelig3MCMC.R b/R/MCMC/zelig3MCMC.R
index 13d4dce..b345df5 100644
--- a/R/MCMC/zelig3MCMC.R
+++ b/R/MCMC/zelig3MCMC.R
@@ -12,14 +12,15 @@
data=eval(out$data))
out$terms <- attr(out$model, "terms")
attr(out$terms,"intercept") <- 0
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
class(out) <- "MCMCZelig"
out
}
zelig3MCMClogit <- zelig3MCMCoprobit <- zelig3MCMCpoisson <-
- zelig3MCMCmnl <- zelig3MCMCregress <-
- zelig3MCMCtobit <- function(res, fcall=NULL, zcall=NULL) {
+ zelig3MCMCmnl <- zelig3MCMCregress <- function(res, fcall=NULL, zcall=NULL) {
out <- list()
out$coefficients <- res
@@ -29,6 +30,9 @@ zelig3MCMClogit <- zelig3MCMCoprobit <- zelig3MCMCpoisson <-
out$model <- model.frame(formula=eval(out$formula),
data=eval(out$data))
out$terms <- attr(out$model, "terms")
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
+
class(out) <- "MCMCZelig"
out
@@ -54,13 +58,38 @@ zelig3MCMCprobit <- function(res, fcall=NULL, zcall=NULL) {
out$model <- model.frame(formula=eval(out$formula),data=eval(out$data))
out$terms <- attr(out$model, "terms")
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
+
class(out) <- "MCMCZelig"
out
}
+ zelig3MCMCtobit <- function(res, fcall=NULL, zcall=NULL) {
+
+ out <- list()
+ out$coefficients <- res
+ out$formula <- zcall$formula
+ if (!is.null(zcall$below)) out$below <- zcall$below
+ else out$below <- 0
+
+ if (!is.null(zcall$above)) out$above <- zcall$above
+ else out$above <- Inf
+
+ out$data <- zcall$data
+
+ out$model <- model.frame(formula=eval(out$formula),
+ data=eval(out$data))
+ out$terms <- attr(out$model, "terms")
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
+ class(out) <- "MCMCZelig"
+
+ out
+}
zelig3MCMCfactanal <- zelig3MCMCordfactanal <- zelig3MCMCmixfactanal <- function(res, fcall=NULL, zcall=NULL) {
@@ -73,6 +102,29 @@ zelig3MCMCprobit <- function(res, fcall=NULL, zcall=NULL) {
data=eval(out$data))
out$terms <- attr(out$model, "terms")
attr(out$terms,"intercept") <- 0
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
+
+ class(out) <- "MCMCZelig"
+
+ out
+}
+
+
+ zelig3MCMCirt1d <- zelig3MCMCirtKd <- function(res, fcall=NULL, zcall=NULL) {
+
+ out <- list()
+ out$coefficients <- res
+ out$formula <- zcall$formula
+ out$data <- zcall$data
+
+ out$model <- model.frame(formula=eval(out$formula),
+ data=eval(out$data))
+ out$terms <- attr(out$model, "terms")
+ attr(out$terms,"intercept") <- 0
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
+
class(out) <- "MCMCZelig"
out
diff --git a/R/class.ind.R b/R/class.ind.R
new file mode 100644
index 0000000..9c11983
--- /dev/null
+++ b/R/class.ind.R
@@ -0,0 +1,11 @@
+#borrow from library(nnet)
+class.ind<-function (cl, levels.data=NULL)
+{
+ n <- length(cl)
+ cl <- as.factor(cl)
+ levels(cl)<-levels.data
+ x <- matrix(0, n, length(levels(cl)))
+ x[(1:n) + n * (unclass(cl) - 1)] <- 1
+ dimnames(x) <- list(names(cl), levels(cl))
+ x
+}
diff --git a/R/gee/param.gee.R b/R/gee/param.gee.R
new file mode 100644
index 0000000..d0f4554
--- /dev/null
+++ b/R/gee/param.gee.R
@@ -0,0 +1,34 @@
+param.gee <- function(object, num = NULL, bootstrap = FALSE) {
+ if (!bootstrap) {
+ coef <- mvrnorm(num, mu=coef(object), Sigma=vcov(object))
+ if (object$zelig == "normal.gee") {
+ sig2 <- object$scale
+ alpha <- sqrt(sig2)
+ res <- cbind(coef, alpha)
+ }
+ else if (object$zelig == "gamma.gee") {
+ alpha<-1/object$scale
+ res <- cbind(coef, alpha)
+ }
+ else
+ res <- coef
+ }
+ else {
+ coef <- coef(object)
+ if (object$family$family == "gaussian") {
+ alpha <- sum(object$residuals^2)/length(object$residuals)
+ res <- c(coef, alpha)
+ }
+ else if (object$family$family == "Gamma") {
+ alpha <- object$scale
+ res <- c(coef, alpha)
+ }
+ else
+ res <- coef
+ }
+ res
+}
+
+
+
+
diff --git a/R/gee/qi.gee.R b/R/gee/qi.gee.R
new file mode 100644
index 0000000..1137f84
--- /dev/null
+++ b/R/gee/qi.gee.R
@@ -0,0 +1,69 @@
+qi.gee <- function(object, simpar, x, x1 = NULL, y = NULL) {
+ check <- FALSE
+ model <- object$zelig
+ k <- length(object$coef)
+ coef <- simpar[,1:k]
+ if (k < ncol(simpar))
+ alpha <- simpar[,(k+1):ncol(simpar)]
+ eta <- coef %*% t(x)
+ theta <- matrix(object$family$linkinv(eta), nrow = nrow(coef))
+ pr <- ev <- matrix(NA, nrow = nrow(theta), ncol = ncol(theta))
+ dimnames(pr) <- dimnames(ev) <- dimnames(theta)
+ if (model =="logit.gee") {
+ check <- TRUE
+ ev <- theta
+ for (i in 1:ncol(theta))
+ pr[,i] <- as.character(rbinom(length(ev[,i]), 1, ev[,i]))
+ }
+ else if (model == "normal.gee") {
+ ev <- theta
+ for (i in 1:nrow(ev))
+ pr[i,] <- rnorm(length(ev[i,]), mean = ev[i,], sd = alpha[i])
+ }
+ else if (model == "gamma.gee") {
+ ev <- theta * 1/alpha
+ for (i in 1:nrow(ev))
+ pr[i,] <- rgamma(length(ev[i,]), shape = theta[i,], scale = 1/alpha[i])
+ }
+ else if (model == "poisson.gee") {
+ ev <- theta
+ for (i in 1:ncol(ev))
+ pr[,i] <- rpois(length(ev[,i]), lambda = ev[,i])
+ }
+ qi <- list(ev = ev, pr = pr)
+ qi.name <- list(ev = "Expected Values: E(Y|X)",
+ pr = "Predicted Values: Y|X")
+ if (!is.null(x1)){
+ theta1 <- matrix(object$family$linkinv(coef %*% t(as.matrix(x1))),
+ nrow = nrow(coef))
+ if (model == "gamma")
+ ev1 <- theta1 * 1/alpha
+ else
+ ev1 <- theta1
+ qi$fd <- ev1-ev
+ qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
+ if (model %in% c("logit", "probit", "relogit")) {
+ qi$rr <- ev1/ev
+ qi.name$rr <- "Risk Ratios: P(Y=1|X1)/P(Y=1|X)"
+ }
+ }
+ if (!is.null(y)) {
+ yvar <- matrix(rep(y, nrow(simpar)), nrow = nrow(simpar), byrow = TRUE)
+ tmp.ev <- yvar - qi$ev
+ if (check)
+ tmp.pr <- yvar - as.integer(qi$pr)
+ else
+ tmp.pr <- yvar - qi$pr
+ qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(simpar))
+ qi$ate.pr <- matrix(apply(tmp.pr, 1, mean), nrow = nrow(simpar))
+ qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
+ qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
+ }
+ list(qi=qi, qi.name=qi.name)
+}
+
+
+
+
+
+
diff --git a/R/gee/vcov.gee.R b/R/gee/vcov.gee.R
new file mode 100644
index 0000000..c0b7543
--- /dev/null
+++ b/R/gee/vcov.gee.R
@@ -0,0 +1,3 @@
+vcov.gee <- function(object, ...) {
+ return(object$robust.variance)
+}
diff --git a/R/gee/zelig2gamma.gee.R b/R/gee/zelig2gamma.gee.R
new file mode 100644
index 0000000..04dec9b
--- /dev/null
+++ b/R/gee/zelig2gamma.gee.R
@@ -0,0 +1,8 @@
+zelig2gamma.gee <- function(formula, model, data, M, ...) {
+ require(gee) || stop("install gee using...")
+ mf <- match.call(expand.dots = TRUE)
+ mf$model <- mf$M <- NULL
+ mf[[1]] <- as.name("gee")
+ mf$family <- as.name("Gamma")
+ as.call(mf)
+}
diff --git a/R/gee/zelig2logit.gee.R b/R/gee/zelig2logit.gee.R
new file mode 100644
index 0000000..0932882
--- /dev/null
+++ b/R/gee/zelig2logit.gee.R
@@ -0,0 +1,8 @@
+zelig2logit.gee <- function(formula, model, data, M, ...) {
+ require(gee) || stop("install gee using...")
+ mf <- match.call(expand.dots = TRUE)
+ mf$model <- mf$M <- NULL
+ mf[[1]] <- as.name("gee")
+ mf$family <- as.name("binomial")
+ as.call(mf)
+}
diff --git a/R/gee/zelig2normal.gee.R b/R/gee/zelig2normal.gee.R
new file mode 100644
index 0000000..c156d78
--- /dev/null
+++ b/R/gee/zelig2normal.gee.R
@@ -0,0 +1,8 @@
+zelig2normal.gee <- function(formula, model, data, M, ...) {
+ require(gee) || stop("install gee using...")
+ mf <- match.call(expand.dots = TRUE)
+ mf$model <- mf$M <- NULL
+ mf[[1]] <- as.name("gee")
+ mf$family <- as.name("gaussian")
+ as.call(mf)
+}
diff --git a/R/gee/zelig2poisson.gee.R b/R/gee/zelig2poisson.gee.R
new file mode 100644
index 0000000..f0cd077
--- /dev/null
+++ b/R/gee/zelig2poisson.gee.R
@@ -0,0 +1,8 @@
+zelig2poisson.gee <- function(formula, model, data, M, ...) {
+ require(gee) || stop("install gee using...")
+ mf <- match.call(expand.dots = TRUE)
+ mf$model <- mf$M <- NULL
+ mf[[1]] <- as.name("gee")
+ mf$family <- as.name("poisson")
+ as.call(mf)
+}
diff --git a/R/load.first.R b/R/load.first.R
index 52a3c9b..e5ae522 100644
--- a/R/load.first.R
+++ b/R/load.first.R
@@ -1,7 +1,11 @@
.onAttach <- function(...) {
- cat("\nPlease refer to http://gking.harvard.edu/zelig for full documentation \n",
- "or help.zelig() for help with commands and models supported by Zelig.\n\n",
- sep='')
+ mylib <- dirname(system.file(package = "Zelig"))
+ ver <- packageDescription("Zelig", lib = mylib)$Version
+ builddate <- packageDescription("Zelig", lib = mylib)$Date
+ cat(paste("## \n## Zelig (Version ", ver, ", built: ", builddate, ")\n", sep = ""))
+ cat("## Please refer to http://gking.harvard.edu/zelig for full documentation \n",
+ "## or help.zelig() for help with commands and models supported by Zelig.\n##\n",
+ sep="")
if(!any(search()=="package:MASS"))
require(MASS)
if(!any(search()=="package:boot"))
diff --git a/R/mcmcei.R b/R/mcmcei.R
new file mode 100644
index 0000000..9b11704
--- /dev/null
+++ b/R/mcmcei.R
@@ -0,0 +1,15 @@
+mcmcei <- function(formula, data, ...) {
+ #if (is.null(rownames(data))) {
+ # rownames(data) <-1:nrow(data)
+ # assign(data, as.character(mc$data), env = .GlobalEnv)
+ #}
+ res <- NULL
+ vars <- model.frame(formula, data)
+ vars <- as.matrix(vars)
+ res$c0 <- vars[,1]
+ res$c1 <- vars[,2]
+ res$r0 <- vars[,3]
+ res$r1 <- vars[,4]
+ res<-as.data.frame(res)
+ res
+}
diff --git a/R/param.MCMCZelig.R b/R/param.MCMCZelig.R
new file mode 100644
index 0000000..007ee15
--- /dev/null
+++ b/R/param.MCMCZelig.R
@@ -0,0 +1,14 @@
+param.MCMCZelig <- function(object, num = NULL, bootstrap = FALSE) {
+ if (bootstrap)
+ stop("For the class of MCMC models, no need to use Bootstrap method.")
+ else
+ res <- object$coefficients
+
+
+ res
+
+}
+
+
+
+
diff --git a/R/param.relogit.R b/R/param.relogit.R
index 44a949f..eaa465c 100644
--- a/R/param.relogit.R
+++ b/R/param.relogit.R
@@ -1,50 +1,22 @@
-param.relogit <- function(object, num, x, bootstrap = FALSE, bootfn = NULL, ...) {
- tau <- eval(object$call$tau, sys.parent())
- if (is.null(object$call$bias.correct))
- object$call$bias.correct <- TRUE
- if (length(tau) == 2) { # without full population information
- tmp0 <- tmp1 <- object
- tmp0$coefficients["(Intercept)"] <- object$correct[1]
- tmp1$coefficients["(Intercept)"] <- object$correct[2]
- pping <- function(object, x, x1, num, bootstrap, bootfn,...) {
- if(!bootstrap) {
- par0 <- param.default(tmp0, num=num, bootstrap=bootstrap)
- par1 <- param.default(tmp1, num=num, bootstrap=bootstrap)
- }
- else {
- dta <- eval(object$data, sys.parent())
- dta <- dta[complete.cases(model.frame(object, dta)),]
- if (is.null(bootfn)) {
- bootfn <- function(data, i, obj) {
- d <- data[i,]
- obj$call$data <- d
- fit <- eval(obj$call, sys.parent())
- fit$coefficients["(Intercept)"] <- fit$correct[1]
- coef0 <- fit$coefficients
- fit$coefficients["(Intercept)"] <- fit$correct[2]
- coef1 <- fit$coefficients
- return(c(coef0, coef1))
- }
- }
- pars <- boot(dta, bootfn, R=num, obj = object, ...)
- colnames(pars$t) <- names(pars$t0)
- par0 <- pars$t[, 1:length(object$coef)]
- par1 <- pars$t[, (length(object$coef)+1):nrow(pars$t)]
- }
- sim00 <- qi.glm(tmp0, par0, x = x)
- P00 <- as.matrix(sim00$qi$ev)
- sim10 <- qi.glm(tmp1, par1, x = x)
- P10 <- as.matrix(sim10$qi$ev)
+param.relogit <- function(object, num, x, bootstrap = FALSE) {
+ if ("relogit2" %in% class(object)) {
+ pping <- function(tmp0, tmp1, num, bootstrap, x) {
+ par0 <- param.relogit(tmp0, num=num, x=x, bootstrap=bootstrap)
+ par1 <- param.relogit(tmp1, num=num, x=x, bootstrap=bootstrap)
+ P00 <- as.matrix(qi.relogit(tmp0, par0, x=x)$qi$ev)
+ P10 <- as.matrix(qi.relogit(tmp1, par1, x=x)$qi$ev)
test <- P00[,1] < P10[,1]
par0 <- as.matrix(par0[test,])
par1 <- as.matrix(par1[test,])
list(par0 = par0, par1 = par1)
}
- tmp <- pping(object, bootstrap=bootstrap, bootfn=bootfn, x=x, x1=x1, num=num, ...)
+ tmp0 <- object$lower.estimate
+ tmp1 <- object$upper.estimate
+ tmp <- pping(tmp0, tmp1, num = num, bootstrap=bootstrap, x=x)
par0 <- tmp$par0
par1 <- tmp$par1
while (nrow(par0) < num) {
- tmp <- pping(object, bootstrap=bootstrap, bootfn=bootfn, x=x, x1=x1, num=num,...)
+ tmp <- pping(tmp0, tmp1, num=num, bootstrap=bootstrap, x=x)
par0 <- rbind(par0, tmp$par0)
par1 <- rbind(par1, tmp$par1)
}
@@ -56,85 +28,11 @@ param.relogit <- function(object, num, x, bootstrap = FALSE, bootfn = NULL, ...)
par1 <- as.matrix(par1)
rownames(par0) <- 1:nrow(par0)
rownames(par1) <- 1:nrow(par1)
- return(list(par0 = par0, par1 = par1))
- }
- else { # with precise population info, or no population info
- object$coefficients["(Intercept)"] <- object$correct
+ return(list(par0 = par0, par1 = par1))
+ } else {
if (!bootstrap)
- simpar <- param.default(object, num=num, bootstrap = bootstrap)
- else {
- tt <- terms(object)
- dta <- eval(object$data, sys.parent())
- dta <- dta[complete.cases(model.frame(tt, dta)),]
- if(is.null(bootfn)) {
- bootfn <- function(data, i, object) {
- d <- data[i,]
- object$call$data <- d
- fit <- eval(object$call, sys.parent())
- fit$coefficients["(Intercept)"] <- correct
- fit$coefficients
- }
- }
- res <- boot(dta, bootfn, R = num, object = object, ...)
- colnames(res$t) <- names(res$t0)
- simpar <- res$t
- }
- return(simpar)
+ return(mvrnorm(num, mu = coef(object), Sigma = vcov(object)))
+ else
+ return(coef(object))
}
}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/R/print.relogit.R b/R/print.relogit.R
index 433e233..1323247 100644
--- a/R/print.relogit.R
+++ b/R/print.relogit.R
@@ -1,28 +1,11 @@
print.relogit <- function(x, digits = max(3, getOption("digits") - 3),
...) {
- if (is.null(x$call$bias.correct))
- x$call$bias.correct <- TRUE
- cat("\nCall: ", deparse(x$call), "\n\n")
- cat("Coefficients")
- if (is.character(co <- x$contrasts))
- cat(" [contrasts: ", apply(cbind(names(co), co),
- 1, paste, collapse = "="), "]")
- cat(":\n")
- print.default(format(x$coefficients, digits = digits),
- print.gap = 2, quote = FALSE)
- if (!is.null(x$call$tau)) {
- cat("\nIntercepts with prior correction for: \n")
- print.default(x$correct, digits = digits,
- print.gap = 1, quote = FALSE)
- }
- cat("\nDegrees of Freedom:", x$df.null, "Total (i.e. Null); ",
- x$df.residual, "Residual\n")
- cat("Null Deviance:\t ", format(signif(x$null.deviance,
- digits)),
- "\nResidual Deviance:", format(signif(x$deviance,
- digits)),
- "\tAIC:", format(signif(x$aic, digits)), "\n")
- if (x$call$bias.correct)
- cat("Rare events bias correction performed.\n")
+ print.glm(x)
+ if (x$prior.correct)
+ cat("\nPrior correction performed with tau =", x$tau, "\n")
+ if (x$weighting)
+ cat("\nWeighting performed with tau =", x$tau, "\n")
+ if (x$bias.correct)
+ cat("Rare events bias correction performed\n")
invisible(x)
}
diff --git a/R/print.relogit2.R b/R/print.relogit2.R
new file mode 100644
index 0000000..e551d3c
--- /dev/null
+++ b/R/print.relogit2.R
@@ -0,0 +1,7 @@
+print.relogit2 <- function(x, digits = max(3, getOption("digits") - 3),
+ ...) {
+ cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
+ print.relogit(x$lower.estimate)
+ print.relogit(x$upper.estimate)
+}
+
diff --git a/R/print.summary.MCMCZelig.R b/R/print.summary.MCMCZelig.R
new file mode 100644
index 0000000..5c0e4c0
--- /dev/null
+++ b/R/print.summary.MCMCZelig.R
@@ -0,0 +1,13 @@
+print.summary.MCMCZelig <- function(x, digits=max(3, getOption("digits") -
+3), ...) {
+ cat("\nCall: ")
+ print(x$call)
+ cat("\n", "Iterations = ", x$start, ":", x$end, "\n", sep = "")
+ cat("Thinning interval =", x$thin, "\n")
+ cat("Number of chains =", x$nchain, "\n")
+ cat("Sample size per chain =", (x$end -
+ x$start)/x$thin + 1, "\n")
+ cat("\n", "Mean, standard deviation, and quantiles for marginal posterior distributions.", "\n")
+ print.matrix(round(x$summary, digits=digits))
+ cat("\n")
+}
diff --git a/R/print.summary.glm.robust.R b/R/print.summary.glm.robust.R
index e26134c..e9d0fe7 100644
--- a/R/print.summary.glm.robust.R
+++ b/R/print.summary.glm.robust.R
@@ -3,69 +3,10 @@ print.summary.glm.robust <-
symbolic.cor = x$symbolic.cor,
signif.stars = getOption("show.signif.stars"), ...)
{
- cat("\nCall:\n")
- cat(paste(deparse(x$call), sep="\n", collapse="\n"), "\n\n", sep="")
- cat("Deviance Residuals: \n")
- if(x$df.residual > 5) {
- x$deviance.resid <- quantile(x$deviance.resid,na.rm=TRUE)
- names(x$deviance.resid) <- c("Min", "1Q", "Median", "3Q", "Max")
- }
- print.default(x$deviance.resid, digits=digits, na = "", print.gap = 2)
-
- if(length(x$aliased) == 0) {
- cat("\nNo Coefficients\n")
- } else {
- ## df component added in 1.8.0
- if (!is.null(df<- x$df) && (nsingular <- df[3] - df[1]))
- cat("\nCoefficients: (", nsingular,
- " not defined because of singularities)\n", sep = "")
- else cat("\nCoefficients:\n")
- coefs <- x$coefficients
- if(!is.null(aliased <- x$aliased) && any(aliased)) {
- cn <- names(aliased)
- coefs <- matrix(NA, length(aliased), 4,
- dimnames=list(cn, colnames(coefs)))
- coefs[!aliased, ] <- x$coefficients
- }
- printCoefmat(coefs, digits=digits, signif.stars=signif.stars,
- na.print="NA", ...)
- }
- ##
- cat("\n(Dispersion parameter for ", x$family$family,
- " family taken to be ", format(x$dispersion), ")\n\n",
- apply(cbind(paste(format.char(c("Null","Residual"),width=8,flag=""),
- "deviance:"),
- format(unlist(x[c("null.deviance","deviance")]),
- digits= max(5, digits+1)), " on",
- format(unlist(x[c("df.null","df.residual")])),
- " degrees of freedom\n"),
- 1, paste, collapse=" "),
- "AIC: ", format(x$aic, digits= max(4, digits+1)),"\n\n",
- "Number of Fisher Scoring iterations: ", x$iter,
- "\n", sep="")
-
- correl <- x$correlation
- if(!is.null(correl)) {
-# looks most sensible not to give NAs for undefined coefficients
-# if(!is.null(aliased) && any(aliased)) {
-# nc <- length(aliased)
-# correl <- matrix(NA, nc, nc, dimnames = list(cn, cn))
-# correl[!aliased, !aliased] <- x$correl
-# }
- p <- NCOL(correl)
- if(p > 1) {
- cat("\nCorrelation of Coefficients:\n")
- if(is.logical(symbolic.cor) && symbolic.cor) {# NULL < 1.7.0 objects
- print(symnum(correl, abbr.col = NULL))
- } else {
- correl <- format(round(correl, 2), nsmall = 2, digits = digits)
- correl[!lower.tri(correl)] <- ""
- print(correl[-1, -p, drop=FALSE], quote = FALSE)
- }
- }
- }
- cat("\nRobust standard errors computed using", x$robust)
- cat("\n")
- invisible(x)
+ class(x) <- "summary.glm"
+ print(x)
+ cat("\nRobust standard errors computed using", x$robust)
+ cat("\n")
+ invisible(x)
}
diff --git a/R/print.summary.relogit.R b/R/print.summary.relogit.R
new file mode 100644
index 0000000..824ad06
--- /dev/null
+++ b/R/print.summary.relogit.R
@@ -0,0 +1,14 @@
+print.summary.relogit <- function(x, digits = max(3, getOption("digits") - 3),
+ ...){
+ class(x) <- "summary.glm"
+ print(x, digits = digits, ...)
+ if (x$prior.correct)
+ cat("\nPrior correction performed with tau =", x$tau, "\n")
+ if (x$weighting)
+ cat("\nWeighting performed with tau =", x$tau, "\n")
+ if (x$bias.correct)
+ cat("Rare events bias correction performed\n")
+ if (!is.null(x$robust))
+ cat("Robust standard errors computed using", x$robust, "\n")
+ invisible(x)
+}
diff --git a/R/print.summary.relogit2.R b/R/print.summary.relogit2.R
new file mode 100644
index 0000000..8772513
--- /dev/null
+++ b/R/print.summary.relogit2.R
@@ -0,0 +1,6 @@
+print.summary.relogit2 <- function(x, digits = max(3, getOption("digits") - 3),
+ ...){
+ cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
+ print(x$lower.estimate)
+ print(x$upper.estimate)
+}
diff --git a/R/qi.MCMCZelig.R b/R/qi.MCMCZelig.R
new file mode 100644
index 0000000..42cedbd
--- /dev/null
+++ b/R/qi.MCMCZelig.R
@@ -0,0 +1,409 @@
+
+qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
+
+ model <- object$zelig
+ qi <- list()
+ check <- FALSE
+
+ if (model %in% c("logit.bayes", "probit.bayes", "oprobit.bayes", "mlogit.bayes"))
+ check <- TRUE
+
+ if (model %in% c("logit.bayes","probit.bayes", "normal.bayes",
+ "poisson.bayes","tobit.bayes")) {
+
+ if (model == "logit.bayes") {
+ coef <- object$coefficients
+ eta <- coef %*% t(x)
+ pr <- ev <- matrix(NA, nrow = nrow(eta), ncol = ncol(eta))
+ dimnames(pr) <- dimnames(ev) <- dimnames(eta)
+ ev <- 1/(1+exp(-eta))
+ for (i in 1:ncol(ev))
+ pr[,i] <- as.character(rbinom(length(ev[,i]), 1, ev[,i]))
+
+ qi$ev <- ev
+ qi$pr <- pr
+ qi.name <- list(ev = "Expected Values: E(Y|X)", pr="Predicted
+ Values: Y|X")
+ }
+ else if (model == "probit.bayes") {
+ coef <- object$coefficients
+ eta <- coef %*% t(x)
+ pr <- ev <- matrix(NA, nrow = nrow(eta), ncol = ncol(eta))
+ dimnames(pr) <- dimnames(ev) <- dimnames(eta)
+ ev <- pnorm(eta)
+ for (i in 1:ncol(ev))
+ pr[,i] <- as.character(rbinom(length(ev[,i]), 1, ev[,i]))
+ qi$ev <- ev
+ qi$pr <- pr
+ qi.name <- list(ev = "Expected Values: E(Y|X)", pr="Predicted
+ Values: Y|X")
+ }
+ else if (model =="normal.bayes") {
+ coef <- object$coefficients[,1:(ncol(object$coefficients)-1)]
+ eta <- coef %*% t(x)
+ ev <- matrix(NA, nrow = nrow(eta), ncol = ncol(eta))
+ dimnames(ev) <- dimnames(eta)
+ ev <- eta
+ qi$ev <- ev
+ qi.name <- list(ev = "Expected Values: E(Y|X)")
+ }
+ else if (model =="tobit.bayes") {
+ coef <- object$coefficients[,1:(ncol(object$coefficients)-1)]
+ sig2 <- object$coefficients[,ncol(object$coefficients)]
+ sig <- sqrt(sig2)
+ eta <- coef %*% t(x)
+ ev <- cev <- matrix(NA, nrow = nrow(eta), ncol = ncol(eta))
+ dimnames(cev) <- dimnames(ev) <- dimnames(eta)
+
+ L2 <- (object$above-eta)/sig
+ L1 <- (object$below-eta)/sig
+ #cev <- eta + sig*(dnorm(L1)-dnorm(L2))/(pnorm(L2)-pnorm(L1))
+ if (object$below==-Inf) temp1<-0
+ else temp1 <- pnorm(L1)*object$below
+ if (object$above==Inf) temp2<-0
+ else temp2 <- (1-pnorm(L2))*object$above
+
+ ev <- temp1+eta*(pnorm(L2)-pnorm(L1))+sig*(dnorm(L1)-dnorm(L2))+temp2
+
+ qi$ev <- ev
+
+ qi.name <- list(ev = "Expected Values: E(Y|X)")
+ }
+ else if (model == "poisson.bayes") {
+ coef <- object$coefficients
+ eta <- coef %*% t(x)
+ pr <- ev <- matrix(NA, nrow = nrow(eta), ncol = ncol(eta))
+ dimnames(pr) <- dimnames(ev) <- dimnames(eta)
+ ev <- exp(eta)
+ for (i in 1:ncol(ev))
+ pr[,i] <- rpois(length(ev[,i]), ev[,i])
+ qi$ev <- ev
+ qi$pr <- pr
+ qi.name <- list(ev = "Expected Values: E(Y|X)", pr="Predicted
+ Values: Y|X")
+
+ }
+
+
+ if (!is.null(x1)) {
+ eta1 <- coef %*% t(x1)
+ if (model == "logit.bayes") {
+ ev1 <- 1/(1+exp(-eta1))
+ rr <-ev1/ev
+ fd <-ev1-ev
+
+ qi$fd <- fd
+ qi$rr <- rr
+
+ qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
+
+ qi.name$rr <- "Risk Ratios: P(Y=1|X1)/P(Y=1|X)"
+ }
+ else if (model == "probit.bayes") {
+ ev1 <- pnorm(eta1)
+ rr <-ev1/ev
+ fd <-ev1-ev
+
+ qi$fd <- fd
+ qi$rr <- rr
+
+ qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
+
+ qi.name$rr <- "Risk Ratios: P(Y=1|X1)/P(Y=1|X)"
+ }
+ else if (model == "normal.bayes") {
+ ev1 <- eta1
+ fd <-ev1-ev
+
+ qi$fd <- fd
+
+ qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
+ }
+ else if (model == "tobit.bayes") {
+ L2 <- (object$above-eta1)/sig
+ L1 <- (object$below-eta1)/sig
+ #cev <- eta + sig*(dnorm(L1)-dnorm(L2))/(pnorm(L2)-pnorm(L1))
+
+ if (object$below==-Inf) temp1<-0
+ else temp1 <- pnorm(L1)*object$below
+ if (object$above==Inf) temp2<-0
+ else temp2 <- (1-pnorm(L2))*object$above
+
+ ev1 <- temp1+eta*(pnorm(L2)-pnorm(L1))+sig*(dnorm(L1)-dnorm(L2))+temp2
+
+ fd <-ev1-ev
+
+ qi$fd <- fd
+
+ qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
+ }
+ else if (model == "poisson.bayes") {
+ ev1 <- exp(eta1)
+ fd <-ev1-ev
+
+ qi$fd <- fd
+
+ qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
+
+ }
+ }
+
+ if (!is.null(y)) {
+ yvar <- matrix(rep(y, nrow(simpar)), nrow = nrow(simpar), byrow = TRUE)
+ tmp.ev <- yvar - qi$ev
+ if (check)
+ tmp.pr <- yvar - as.integer(qi$pr)
+ else
+ tmp.pr <- yvar - qi$pr
+ qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(simpar))
+ qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
+ if (model %in% c("logit", "probit", "poisson"))
+ {
+ qi$ate.pr <- matrix(apply(tmp.pr, 1, mean), nrow = nrow(simpar))
+ qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
+ }
+ }
+
+ list(qi=qi, qi.name=qi.name)
+ }
+ else if ((model =="oprobit.bayes") || (model == "mlogit.bayes")) {
+ if (model == "oprobit.bayes") {
+ library(stats)
+ p <- dim(model.matrix(object, data=eval(object$data)))[2]
+ coef <- object$coefficients[,1:p]
+ eta <- coef %*% t(x)
+
+ level <- ncol(object$coefficients)-p+2
+ gamma<-matrix(NA, nrow(object$coefficients), level+1)
+
+ gamma[,1] <- rep(-Inf, nrow(gamma))
+ gamma[,2] <- rep(0, nrow(gamma))
+ gamma[,ncol(gamma)]<-rep(Inf, nrow(gamma))
+
+ if (ncol(gamma)>3)
+ gamma[,3:(ncol(gamma)-1)] <- object$coefficients[,(p+1):ncol(object$coefficients)]
+
+
+ ev <- array(NA, c(nrow(eta), level, ncol(eta)))
+ pr <- matrix(NA, nrow(eta), ncol(eta))
+# dimnames(pr)[1] <- dimnames(ev)[1] <- dimnames(eta)[1]
+# dimnames(pr)[2] <- dimnames(ev)[3] <- dimnames(eta)[2]
+
+ for (j in 1:level)
+ ev[,j,] <- pnorm(gamma[,j+1]-eta) - pnorm(gamma[,j]-eta)
+
+ colnames(ev) <- levels(model.response(model.frame(object)))
+
+ for (j in 1:nrow(pr)) {
+ mu <- eta[j,]
+# pr[j,]<-as.character(cut(mu, gamma[j,],
+# labels=as.factor(1:level)))
+ pr[j,]<-as.character(cut(mu, gamma[j,], labels=colnames(ev)))
+ }
+
+
+ colnames(ev) <- levels(model.response(model.frame(object)))
+ qi$ev <- ev
+ qi$pr <- pr
+ qi.name <- list(ev = "Expected Values: P(Y=j|X)", pr="Predicted
+ Values: Y|X")
+ }
+ else if (model == "mlogit.bayes") {
+ library(stats)
+ resp <- model.response(model.frame(object))
+ level <- length(table(resp))
+
+ p <- dim(model.matrix(eval(object),data=eval(object$data)))[2]
+
+ coef <- object$coefficients
+
+
+ eta <- array(NA, c(nrow(coef),level, nrow(x)))
+
+ eta[,1,]<-matrix(0, dim(eta)[1],dim(eta)[3])
+ for (j in 2:level) {
+ ind <- (1:p)*(level-1)-(level-j)
+ eta[,j,]<- coef[,ind]%*%t(x)
+ }
+
+ eta<-exp(eta)
+
+ ev <- array(NA, c(nrow(coef), level, nrow(x)))
+ pr <- matrix(NA, nrow(coef), nrow(x))
+
+
+ colnames(ev) <- rep(NA, level)
+
+ for (k in 1:nrow(x))
+ for (j in 1:level)
+ {
+ ev[,j,k] <- eta[,j,k]/rowSums(eta[,,k])
+ }
+
+ for (j in 1:level)
+ {
+ colnames(ev)[j] <- paste("P(Y=", j, ")", sep="")
+ }
+
+ for (k in 1:nrow(x)) {
+ probs <- as.matrix(ev[,,k])
+ temp <- apply(probs, 1, FUN=rmultinom, n=1, size=1)
+ temp <- as.matrix(t(temp)%*%(1:nrow(temp)))
+ pr <- apply(temp,2,as.character)
+ }
+
+ qi$ev <- ev
+ qi$pr <- pr
+ qi.name <- list(ev = "Expected Values: P(Y=j|X)", pr="Predicted
+ Values: Y|X")
+ }
+
+
+ if (!is.null(x1)) {
+
+ if (model == "oprobit.bayes") {
+
+ eta1 <- coef %*% t(x1)
+
+ ev1 <- array(NA, c(nrow(eta), level, ncol(eta)))
+
+ for (j in 1:level)
+ ev1[,j,] <- pnorm(gamma[,j+1]-eta1) - pnorm(gamma[,j]-eta1)
+
+ rr <-ev1/ev
+ fd <-ev1-ev
+
+ qi$fd <- fd
+ qi$rr <- rr
+
+
+ qi.name$fd <- "First Differences in Expected Values: P(Y=j|X1)-P(Y=j|X)"
+
+ qi.name$rr <- "Risk Ratios: P(Y=j|X1)/P(Y=j|X)"
+
+ }
+ else if (model == "mlogit.bayes") {
+
+ eta1 <- array(NA, c(nrow(coef),level, nrow(x1)))
+
+ eta1[,1,]<-matrix(0, dim(eta1)[1],dim(eta1)[3])
+ for (j in 2:level) {
+ ind <- (1:p)*(level-1)-(level-j)
+ eta1[,j,]<- coef[,ind]%*%t(x1)
+ }
+
+ eta1<-exp(eta1)
+
+ ev1 <- array(NA, c(nrow(eta1), level, nrow(x1)))
+
+
+ for (k in 1:nrow(x))
+ for (j in 1:level)
+ {
+ ev1[,j,k] <- eta1[,j,k]/rowSums(eta1[,,k])
+
+ }
+
+ rr <-ev1/ev
+ fd <-ev1-ev
+
+ qi$fd <- fd
+ qi$rr <- rr
+
+
+ qi.name$fd <- "First Differences in Expected Values: P(Y=j|X1)-P(Y=j|X)"
+
+ qi.name$rr <- "Risk Ratios: P(Y=j|X1)/P(Y=j|X)"
+
+ }
+ }
+
+ if (!is.null(y)) {
+ yvar <- matrix(rep(y, nrow(simpar)), nrow = nrow(simpar), byrow =
+ TRUE)
+
+ levels.names<-levels(as.factor(y))
+ yvar1<- pr1 <-array(NA, c(nrow(yvar), level, ncol(yvar)))
+
+
+ for (j in 1:nrow(yvar))
+ {
+ yvar1[j,,]<-t(class.ind(yvar[j,], levels.names))
+ if (check)
+ pr1[j,,]<-t(class.ind(as.integer(qi$pr[j,]), levels.names))
+ else
+ pr1[j,,]<-t(class.ind(qi$pr[j,],levels.names))
+ }
+
+ tmp.ev <- yvar1 - qi$ev
+ tmp.pr <- yvar1 - pr1
+ qi$ate.ev <- matrix(apply(tmp.ev, 2, rowMeans), nrow = nrow(simpar))
+ qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
+ if (model %in% c("oprobit.bayes", "mlogit.bayes"))
+ {
+ qi$ate.pr <- matrix(apply(tmp.pr, 2, rowMeans), nrow = nrow(simpar))
+ qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
+ }
+ }
+
+ list(qi=qi, qi.name=qi.name)
+
+
+
+ }
+ else if (model == "ei.hier" || model == "ei.dynamic") {
+ if (!any(class(x)=="cond")) stop("set 'cond=TRUE' in setx.\n")
+ else
+ {
+ coef <- object$coefficients
+ n <- nrow(x)
+ if (is.null(object$N))
+ N<-rep(1,n)
+ else N <- eval(object$N)
+ ev <- array(NA, c(nrow = nrow(coef), 2,2, n))
+ pr <- array(NA, c(nrow = nrow(coef), 2,2, n))
+ nlen<-length(coef[,1])
+ for (j in 1:2) {
+ ev[,j,1,] <- t(apply(coef[,((1:n)+(j-1)*n)],
+ 1,"*", x[,j])*N)
+ ev[,j,2,] <- t(apply((1-coef[,((1:n)+(j-1)*n)]), 1,"*",
+ x[,j])*N)
+ for (i in 1:n)
+ {
+ size<-round(x[i,j]*N[i])
+ pr[,j,1,i] <-rbinom(prob=coef[,(i+(j-1)*n)], n=nlen, size=size)
+
+ pr[,j,2,i] <- x[i,j]*N[i]-pr[,j,1,i]
+ }
+
+ }
+# dimnames(ev)[[1]] <- dimnames(pr)[[4]] <- 1:nrow(coef)
+ dimnames(ev)[[4]] <- dimnames(pr)[[4]] <- rownames(x)
+ dimnames(ev)[[2]] <- dimnames(pr)[[2]] <- colnames(x)
+ dimnames(ev)[[3]] <- dimnames(pr)[[3]] <- colnames(model.response(object$model))
+ class(ev) <- class(pr) <- c("ei", "array")
+ qi$ev <- ev
+ qi$pr <- pr
+ qi.name <- list(ev = "Expected In sample predictions at aggregate
+ level", pr = "In sample predictions at aggregate level")
+
+ }
+
+ list(qi=qi, qi.name=qi.name)
+ }
+ else if ( model %in% c("factor.bayes", "factor.ord",
+ "factor.mix", "irt1d", "irtkd"))
+ {
+ stop("sim procedure not applicable since no explanatory
+ variables are involved.\n")
+ list(qi=qi)
+ }
+
+
+
+}
+
+
+
+
+
diff --git a/R/qi.relogit.R b/R/qi.relogit.R
index 59bdb40..9cc4324 100644
--- a/R/qi.relogit.R
+++ b/R/qi.relogit.R
@@ -1,14 +1,11 @@
qi.relogit <- function(object, simpar, x, x1 = NULL, y = NULL) {
- tau <- eval(object$call$tau, sys.parent())
- if (is.null(object$call$bias.correct))
- object$call$bias.correct <- TRUE
- num <- nrow(simpar$par0)
- tmp0 <- tmp1 <- object
- tmp0$coefficients["(Intercept)"] <- object$correct[1]
- tmp1$coefficients["(Intercept)"] <- object$correct[2]
- if (length(tau) == 2){
- low <- qi.glm(tmp0, simpar$par0, x, x1)
- up <- qi.glm(tmp1, simpar$par1, x, x1)
+ if ("relogit2" %in% class(object)) {
+ num <- nrow(simpar$par0)
+ tmp0 <- object$lower.estimate
+ tmp1 <- object$upper.estimate
+
+ low <- qi.relogit(tmp0, simpar$par0, x, x1)
+ up <- qi.relogit(tmp1, simpar$par1, x, x1)
PP <- PR <- array(NA, dim = c(num, 2, nrow(x)),
dimnames = list(NULL, c("Lower Bound", "Upper Bound"),
rownames(x)))
@@ -22,8 +19,8 @@ qi.relogit <- function(object, simpar, x, x1 = NULL, y = NULL) {
rownames(x)))
sim01 <- qi.glm(tmp0, par0, x = x1, x1 = NULL)
sim11 <- qi.glm(tmp1, par1, x = x1, x1 = NULL)
- tau0 <- object$call$tau[1]
- tau1 <- object$call$tau[2]
+ tau0 <- object$lower.estimate$tau
+ tau1 <- object$upper.estimate$tau
P01 <- as.matrix(sim01$qi$ev)
P11 <- as.matrix(sim11$qi$ev)
OR <- (P10/(1-P10)) / (P00/(1-P00))
@@ -34,8 +31,8 @@ qi.relogit <- function(object, simpar, x, x1 = NULL, y = NULL) {
RD <- as.matrix((sqrt(OR)-1) / (sqrt(OR)+1))
## checking monotonicity
y.bar <- mean(object$y)
- beta0.e <- tmp0$coefficients
- beta1.e <- tmp1$coefficients
+ beta0.e <- coef(tmp0)
+ beta1.e <- coef(tmp1)
## evaluating RD at tau0 and tau1
RD0.p <- 1/(1+exp(-t(beta0.e) %*% t(x1))) - 1/(1+exp(-t(beta0.e) %*% t(x)))
RD1.p <- 1/(1+exp(-t(beta1.e) %*% t(x1))) - 1/(1+exp(-t(beta1.e) %*% t(x)))
@@ -71,9 +68,8 @@ qi.relogit <- function(object, simpar, x, x1 = NULL, y = NULL) {
qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
}
- out <- list(qi = qi, qi.name = qi.name)
+ return(list(qi = qi, qi.name = qi.name))
}
- else
- out <- qi.glm(object = object, simpar = simpar, x = x, x1 = x1, y = y)
- out
+ else
+ return(qi.glm(object = object, simpar = simpar, x = x, x1 = x1, y = y))
}
diff --git a/R/qi.survreg.R b/R/qi.survreg.R
index fea4e78..f9c78c2 100644
--- a/R/qi.survreg.R
+++ b/R/qi.survreg.R
@@ -1,71 +1,86 @@
-qi.survreg <- function(object, simpar, x, x1 = NULL, y = NULL, cond.data = NULL) {
+qi.survreg <- function(object, simpar, x, x1 = NULL, y = NULL) {
model <- object$zelig
k <- length(object$coef)
sim.coef <- as.matrix(simpar[,1:k])
- if (model %in% c("weibull", "Weibull", "lognorm")) #{
+ if (model %in% c("weibull", "Weibull")) {
+ if (ncol(simpar) == (k + 1))
+ sim.scale <- simpar[,(k+1):ncol(simpar)]
+ else
+ sim.scale <- rep(object$scale, nrow(simpar))
+ }
+ else if (model == "lognorm")
sim.scale <- simpar[,(k+1):ncol(simpar)]
if (!is.null(y)) {
- status <- x[,1]
- x <- as.data.frame(x[,2:ncol(x)])
+ status <- y[,2]
+ y <- y[,1]
}
link <- survreg.distributions[[object$dist]]$itrans
- ev.pr.surv <- function(sim.coef, sim.scale, x, link) {
+ ev.surv <- function(model, sim.coef, sim.scale, x, link) {
eta <- sim.coef %*% t(x)
theta <- as.matrix(apply(eta, 2, link))
- ev <- pr <- matrix(NA, ncol=ncol(theta), nrow=nrow(theta))
- dimnames(pr) <- dimnames(ev) <- dimnames(theta)
- if (model == "exp") {
- ev <- theta
- for (i in 1:nrow(pr))
- pr[i,] <- rexp(length(ev[i,]), rate = 1/ev[i,])
+ if (model == "lognorm") {
+ ev <- exp(log(theta) + 0.5*(exp(sim.scale))^2)
+ dimnames(ev) <- dimnames(theta)
}
else if (model == "weibull" || model == "Weibull") {
ev <- theta * gamma(1 + exp(sim.scale))
- for (i in 1:nrow(pr))
- pr[i,] <- rweibull(length(ev[i,]), shape=1/exp(sim.scale[i]),
- scale=theta[i,])
+ dimnames(ev) <- dimnames(theta)
}
- else if (model == "lognorm") {
- ev <- exp(log(theta) + 0.5*(1/exp(sim.scale))^2)
- for (i in 1:nrow(pr))
- pr[i,] <- rlnorm(length(ev[i,]), meanlog = log(theta[i,]),
- sdlog = 1/exp(sim.scale[i]))
- }
- list(ev=ev, pr=pr)
- }
- qi <- ev.pr.surv(sim.coef, sim.scale, x, link)
- if (!is.null(y)) {
- if (is.null(cond.data))
- stop("`cond.data' is required for the exponential, Weibull, and lognormal models.")
- call <- object$call
- call$by <- NULL
- call$data <- cond.data
- est0 <- eval(call, sys.frame())
- par0 <- matrix(param(est0, bootstrap = TRUE), nrow = 1)
- coef0 <- matrix(par0[,1:k], nrow = 1)
- if (model %in% c("weibull", "Weibull", "lognorm"))
- scale0 <- rep(par0[,(k+1):ncol(par0)], nrow(x))
- qi0 <- ev.pr.surv(coef0, scale0, x, link)
- y1 <- y
- tmp <- status
- y[status == 0] <- qi0$pr[(status == 0), 1]
- status[status == 0] <- as.integer(y1[status == 0] < y[status == 0])
- while (sum(status) < length(status)) {
- qi0 <- ev.pr.surv(coef0, scale0, x, link)
- y[status == 0] <- qi0$pr[status == 0, 1]
- status[status == 0] <- as.integer(y1[status == 0] < y[status == 0])
+ else if (model == "exp") {
+ ev <- theta
}
+ list(ev = as.matrix(ev), theta = as.matrix(theta))
+ }
+ pr.surv <- function(model, theta, sim.scale, ev) {
+ if (model == "exp")
+ pr <- rexp(length(ev), rate = 1/ev)
+ else if (model == "weibull" || model == "Weibull")
+ pr <- rweibull(length(ev), shape=1/exp(sim.scale),
+ scale=theta)
+ else if (model == "lognorm")
+ pr <- rlnorm(length(ev), meanlog = log(theta),
+ sdlog = exp(sim.scale))
+ pr
}
+ ev <- ev.surv(model, sim.coef, sim.scale, x, link)
+ pr <- matrix(NA, ncol=ncol(ev$ev), nrow=nrow(ev$ev))
+ dimnames(pr) <- dimnames(ev$ev)
+ for (i in 1:nrow(ev$ev))
+ pr[i,] <- pr.surv(model, ev$theta[i,], sim.scale[i], ev$ev[i,])
+ qi <- list(ev = ev$ev, pr = pr)
qi.name <- list(ev="Expected Values: E(Y|X)",
pr="Predicted Values: Y|X")
if (!is.null(x1)) {
- eta1.sim <- sim.coef %*% t(x1)
- ev1 <- apply(eta1.sim, 2, link)
- qi$fd <- ev1 - qi$ev
+ ev1 <- ev.surv(model, sim.coef, sim.scale, x1, link)
+ qi$fd <- ev1$ev - qi$ev
qi.name$fd <- "First Differences: E(Y|X1)-E(Y|X)"
}
if (!is.null(y)) {
- yvar <- matrix(rep(y, nrow(simpar)), nrow = nrow(simpar), byrow = TRUE)
+ tmp <- list(ev = ev$ev[, which(status == 0)], theta = ev$theta[, which(status == 0)])
+ y.obs <- matrix(y[status == 1], nrow = nrow(qi$ev),
+ ncol = length(y[status == 1]), byrow = TRUE)
+ y.imp <- matrix(NA, nrow = nrow(qi$ev), ncol = length(y[status == 0]))
+ tmp.scale <- c(matrix(sim.scale, nrow = length(sim.scale),
+ ncol = length(y[status == 0])))
+ y.imp <- matrix(pr.surv(model, tmp$theta, tmp.scale, tmp$ev),
+ nrow = nrow(qi$ev), ncol = length(y[status == 0]))
+ y.c <- y[status == 0]
+ idx <- t(apply(y.imp, 1, '>=', y.c))
+ count <- 1
+ while ((sum(idx) < length(idx)) & count < 1001) {
+ count <- count + 1
+ tmp.idx <- which(!idx, arr.ind = TRUE)
+ y.imp[tmp.idx] <- pr.surv(model, tmp$theta[tmp.idx],
+ sim.scale[tmp.idx[,1]], tmp$ev[tmp.idx])
+ idx[tmp.idx] <- y.imp[tmp.idx] >= y.c[tmp.idx[,2]]
+ }
+ if (count == 1001) {
+ warning(" Maximum number of imputed values (1000) reached for censored Y. \n Using censoring point as observed value, since Pr(Y > Yc | sims) <= 0.001.")
+ y.imp[which(idx == 0, arr.ind = TRUE)] <- y.c[which(idx == 0, arr.ind == TRUE)[,2]]
+ }
+ yvar <- matrix(NA, ncol = length(y), nrow = nrow(qi$ev))
+ yvar[, which(status == 1)] <- y.obs
+ yvar[, which(status == 0)] <- y.imp
tmp.ev <- yvar - qi$ev
tmp.pr <- yvar - qi$pr
qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(simpar))
diff --git a/R/relogit.R b/R/relogit.R
index 24fcb98..6729282 100644
--- a/R/relogit.R
+++ b/R/relogit.R
@@ -1,47 +1,92 @@
-relogit <- function(formula, data=sys.parent(), tau=NULL,
- bias.correct=TRUE, ...){
+relogit <- function(formula, data = sys.parent(), tau = NULL,
+ bias.correct = TRUE, case.control = "prior", ...){
mf <- match.call()
- mf$tau <- mf$bias.correct <- NULL
- if (!is.null(tau))
+ mf$tau <- mf$bias.correct <- mf$case.control <- NULL
+ if (!is.null(tau)) {
tau <- unique(tau)
+ if (length(case.control) > 1)
+ stop("You can only choose one option for case control correction.")
+ ck1 <- grep("p", case.control)
+ ck2 <- grep("w", case.control)
+ if (length(ck1) == 0 & length(ck2) == 0)
+ stop("choose either case.control = \"prior\" or case.control = \"weighting\"")
+ if (length(ck2) == 0)
+ weighting <- FALSE
+ else
+ weighting <- TRUE
+ }
if (length(tau) > 2)
stop("tau must be a vector of length less than or equal to 2")
- mf[[1]] <- glm
- mf$family <- binomial(link="logit")
- res <- eval(as.call(mf))
- res$call <- match.call(expand.dots = TRUE)
- ## prior correction
- y <- res$y
- ybar <- mean(y)
- #w1 <- tau/ybar
- #w0 <- (1-tau)/(1-ybar)
- #wi <- w1*y + w0*(1-y)
- X <- model.matrix(res)
- if (bias.correct){
- pihat <- fitted(res)
- # if (is.null(tau))
- # W <- pihat * (1 - pihat) * wi
- # else
- W <- pihat * (1 - pihat)
- Qdiag <- lm.influence(glm(y ~ X-1, weights=W))$hat
- #if (is.null(tau))
- # xi <- 0.5 * Qdiag * ((1+w0)*pihat-w0)
- #else
- xi <- 0.5 * Qdiag * (2*pihat - 1)
- res$coefficients <- res$coefficients -
- glm(xi ~ X - 1, weights=W)$coefficients
+ else if (length(tau)==2) {
+ mf[[1]] <- relogit
+ res <- list()
+ mf$tau <- min(tau)
+ res$lower.estimate <- eval(as.call(mf))
+ mf$tau <- max(tau)
+ res$upper.estimate <- eval(as.call(mf))
+ class(res) <- c("relogit2", "relogit")
+ return(res)
}
- if (!is.null(tau)){
- if (tau <= 0 || tau >= 1)
- stop("\ntau needs to be between 0 and 1.\n")
- res$correct <- res$coefficients["(Intercept)"] -
- log(((1-tau)/tau) * (ybar/(1-ybar)))
- names(res$correct) <- paste("tau =", round(tau, 4))
+ else {
+ mf[[1]] <- glm
+ mf$family <- binomial(link="logit")
+ y2 <- model.response(model.frame(mf$formula, data))
+ if (is.matrix(y2))
+ y <- y2[,1]
+ else
+ y <- y2
+ ybar <- mean(y)
+ if (weighting) {
+ w1 <- tau/ybar
+ w0 <- (1-tau)/(1-ybar)
+ wi <- w1*y + w0*(1-y)
+ mf$weights <- wi
+ }
+ res <- eval(as.call(mf))
+ res$call <- match.call(expand.dots = TRUE)
+ res$tau <- tau
+ X <- model.matrix(res)
+ ## bias correction
+ if (bias.correct){
+ pihat <- fitted(res)
+ if (is.null(tau)) # w_i = 1
+ wi <- rep(1, length(y))
+ else if (weighting)
+ res$weighting <- TRUE
+ else {
+ w1 <- tau/ybar
+ w0 <- (1-tau)/(1-ybar)
+ wi <- w1*y + w0*(1-y)
+ res$weighting <- FALSE
+ }
+ W <- pihat * (1 - pihat) * wi
+ ##Qdiag <- diag(X%*%solve(t(X)%*%diag(W)%*%X)%*%t(X))
+ Qdiag <- lm.influence(lm(y ~ X-1, weights=W))$hat/W
+ if (is.null(tau)) # w_1=1 since tau=ybar
+ xi <- 0.5 * Qdiag * (2*pihat - 1)
+ else
+ xi <- 0.5 * Qdiag * ((1+w0)*pihat-w0)
+ res$coefficients <- res$coefficients -
+ lm(xi ~ X - 1, weights=W)$coefficients
+ res$bias.correct <- TRUE
+ }
+ else
+ res$bias.correct <- FALSE
+ ## prior correction
+ if (!is.null(tau) & !weighting){
+ if (tau <= 0 || tau >= 1)
+ stop("\ntau needs to be between 0 and 1.\n")
+ res$coefficients["(Intercept)"] <- res$coefficients["(Intercept)"] -
+ log(((1-tau)/tau) * (ybar/(1-ybar)))
+ res$prior.correct <- TRUE
+ }
+ else
+ res$prior.correct <- FALSE
+ res$linear.predictors <- t(res$coefficients) %*% t(X)
+ res$fitted.values <- 1/(1+exp(-res$linear.predictors))
+ res$zelig <- "relogit"
+ class(res) <- c("relogit", "glm")
+ return(res)
}
- res$linear.predictors <- t(res$coef) %*% t(X)
- res$fitted.values <- 1/(1+exp(-res$linear.predictors))
- res$zelig <- "relogit"
- class(res) <- "relogit"
- res
}
diff --git a/R/setx.relogit2.R b/R/setx.relogit2.R
new file mode 100644
index 0000000..b554175
--- /dev/null
+++ b/R/setx.relogit2.R
@@ -0,0 +1,7 @@
+setx.relogit2 <- function(object, fn = list(numeric = mean,
+ ordered = median, other = mode), data =
+ NULL, cond = FALSE,
+ counter = NULL, ...) {
+ return(setx.default(object$lower.estimate))
+}
+
diff --git a/R/sim.cond.R b/R/sim.cond.R
index 86862cf..a46dc25 100644
--- a/R/sim.cond.R
+++ b/R/sim.cond.R
@@ -1,7 +1,6 @@
sim.cond <- function(object, x, x1=NULL, num=c(1000, 100),
qoi = c("ev", "pr"), prev = NULL,
- bootstrap = FALSE, bootfn=NULL,
- cond.data = NULL, ...) {
+ bootstrap = FALSE, bootfn=NULL, ...) {
if (!is.null(x1)) {
warning("First Differences are not calculated in conditional prediction models.")
x1 <- NULL
@@ -40,15 +39,11 @@ sim.cond <- function(object, x, x1=NULL, num=c(1000, 100),
else
simpar <- prev
}
- if (class(object)[1] == "survreg")
- simqi <- qi(object, simpar = simpar, x = xvar, x1 = x1, y = yvar,
- cond.data = cond.data)
- else
- simqi <- qi(object, simpar = simpar, x = xvar, x1 = x1, y = yvar)
+ simqi <- qi(object, simpar = simpar, x = xvar, x1 = x1, y = yvar)
class(xvar) <- c("matrix", "cond")
- c <- match.call()
- c$num <- num
- res <- list(x=xvar, x1=x1, call = c, zelig.call = object$call,
+ ca <- match.call()
+ ca$num <- num
+ res <- list(x=xvar, x1=x1, call = ca, zelig.call = object$call,
par = simpar, qi=simqi$qi, qi.name=simqi$qi.name)
class(res) <- "zelig"
res
diff --git a/R/sim.default.R b/R/sim.default.R
index 876621b..f81c0b6 100644
--- a/R/sim.default.R
+++ b/R/sim.default.R
@@ -13,10 +13,10 @@ sim.default <- function(object, x, x1=NULL, num=c(1000, 100),
num <- num[2]
}
if (is.null(prev)) {
- if (!bootstrap & any(class(object) != "relogit"))
+ if (any(class(object) == "relogit"))
+ simpar <- param.relogit(object, num=num, x=x, bootstrap=bootstrap)
+ else if (!bootstrap)
simpar <- param(object, num=num, bootstrap=bootstrap)
- else if (any(class(object) == "relogit"))
- simpar <- param.relogit(object, num=num, x=x, bootstrap=bootstrap, bootfn=bootfn, ...)
else {
tt <- terms(object)
dta <- eval(object$data, sys.parent())
diff --git a/R/summarize.ei.R b/R/summarize.ei.R
new file mode 100644
index 0000000..2436c67
--- /dev/null
+++ b/R/summarize.ei.R
@@ -0,0 +1,24 @@
+summarize.ei <- function(x, rows = NULL, cip, stats, subset = NULL) {
+ if (is.function(subset)) { # subset = all; all is class "function"
+ res <- apply(x, c(2,3,4), summarize.default, stats = stats, cip = cip)
+ dimnames(res)[[4]] <- rows
+ }
+ if (is.null(subset)){# subset = NULL; summarizes all obs at once
+ tmp <- NULL
+ tmp <- apply(x, c(2,3), rbind, tmp)
+ res <- apply(tmp, c(2,3), summarize.default,
+ stats = stats, cip = cip)
+ }
+ if (is.numeric(subset)) { # subset=integer, summarizes identified obs
+ if (length(subset) > 1) {
+ res <- apply(x[, , , subset], c(2,3,4), summarize.default,
+ stats = stats, cip = cip)
+ dimnames(res)[[4]] <- rows
+ }
+ else
+ res <- apply(x[, , subset], 2, summarize.default,
+ stats = stats, cip = cip)
+ }
+ dimnames(res)[2:3] <- dimnames(x)[2:3]
+ res
+}
diff --git a/R/summary.MCMCZelig.R b/R/summary.MCMCZelig.R
new file mode 100644
index 0000000..c297013
--- /dev/null
+++ b/R/summary.MCMCZelig.R
@@ -0,0 +1,19 @@
+
+summary.MCMCZelig <- function(object, quantiles = c(0.025, 0.5, 0.975), ...) {
+ require(coda)
+ out <- list()
+ out$summary <- cbind(summary(object$coefficients)$statistics[,1:2],
+ summary(object$coefficients,
+ quantiles=quantiles)$quantiles)
+
+ colnames(out$summary) <- c("Mean", "SD", paste(quantiles*100, "%",sep=""))
+ stuff <- attributes(object$coefficients)
+ out$call <- object$call
+ out$start <- stuff$mcpar[1]
+ out$end <- stuff$mcpar[2]
+ out$thin <- stuff$mcpar[3]
+ out$nchain <- 1
+ class(out) <- "summary.MCMCZelig"
+ out
+}
+
diff --git a/R/summary.relogit.R b/R/summary.relogit.R
index 32deddb..bce1f7b 100644
--- a/R/summary.relogit.R
+++ b/R/summary.relogit.R
@@ -1,31 +1,24 @@
summary.relogit <- function(object, ...) {
- if (!is.numeric(object$call$tau))
- tau <- eval(object$call$tau, sys.frame(sys.parent()))
- else
- tau <- object$call$tau
- if (is.null(object$call$bias.correct))
- object$call$bias.correct <- TRUE
- dta <- eval(object$call$data, sys.parent())
- tt <- terms(object)
- dta <- dta[complete.cases(model.frame(tt, dta)),]
- n <- nrow(dta)
- k <- ncol(dta)
- bias.cor <- function(res.obj, bc, n, k){
- if (bc) {
- res.obj$cov.unscaled <- res.obj$cov.unscaled * (n/(n+k))^2
- res.obj$cov.scaled <- res.obj$cov.unscaled * res.obj$dispersion
- res.obj$coef[,2] <- sqrt(diag(res.obj$cov.scaled))
- res.obj$coef[,3] <- res.obj$coef[,1] / res.obj$coef[,2]
- res.obj$coef[,4 ] <- 2*pt(-abs(res.obj$coef[,3]), res.obj$df.residual)
- }
- res.obj
+
+ dta <- model.matrix(terms(object), data=model.frame(object))
+ class(object) <- class(object)[2]
+ res <- summary(object, ...)
+ if (object$bias.correct) {
+ n <- nrow(dta)
+ k <- ncol(dta)
+ res$cov.unscaled <- res$cov.unscaled * (n/(n+k))^2
+ res$cov.scaled <- res$cov.unscaled * res$dispersion
+ res$coef[,2] <- sqrt(diag(res$cov.scaled))
+ res$coef[,3] <- res$coef[,1] / res$coef[,2]
+ res$coef[,4 ] <- 2*pt(-abs(res$coef[,3]), res$df.residual)
}
- res <- summary.glm(object, ...)
- res <- bias.cor(res, object$call$bias.correct, n, k)
res$call <- object$call
- res$correct <- object$correct
- class(res) <- "relogit"
- res
+ res$tau <- object$tau
+ res$bias.correct <- object$bias.correct
+ res$prior.correct <- object$prior.correct
+ res$weighting <- object$weighting
+ class(res) <- "summary.relogit"
+ return(res)
}
diff --git a/R/summary.relogit2.R b/R/summary.relogit2.R
new file mode 100644
index 0000000..c6ee09b
--- /dev/null
+++ b/R/summary.relogit2.R
@@ -0,0 +1,21 @@
+summary.relogit2 <- function(object, ...) {
+
+ res <- list()
+ res$lower.estimate <- summary.relogit(object$lower.estimate)
+ res$upper.estimate <- summary.relogit(object$upper.estimate)
+ res$call <- object$call
+ class(res) <- "summary.relogit2"
+ return(res)
+}
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/R/systemfit/callsystemfit.R b/R/systemfit/callsystemfit.R
new file mode 100644
index 0000000..55a592f
--- /dev/null
+++ b/R/systemfit/callsystemfit.R
@@ -0,0 +1,8 @@
+callsystemfit<-function(formula,data,eqns,method,omit=NULL,constrain=NULL,...){
+
+t<-terms.multiple(formula,omit,constrain)
+out<-systemfit(data=data,eqns=eqns,method=method,...)
+out$terms<-t
+class(out)<-c(class(out),"multiple")
+return (out)
+}
diff --git a/R/systemfit/cmsystemfit.R b/R/systemfit/cmsystemfit.R
new file mode 100644
index 0000000..09b5de3
--- /dev/null
+++ b/R/systemfit/cmsystemfit.R
@@ -0,0 +1,19 @@
+cmsystemfit<-function(formu,omit=NULL,...){
+ tr<-terms.multiple(formu,omit)
+ ev<-attr(tr,"term.labels")
+ dv<-all.vars(attr(tr,"variables")[[2]],unique=FALSE)
+ om<-attr(tr,"omit")
+ syst<-res<-list()
+ for(i in 1:length(dv)){
+ syst[[i]]<- ev[om[i,]==0]
+ res[[i]]<-paste(dv[i],"~")
+ for(j in 1:(length(syst[[i]])-1)){
+ res[[i]]<-paste(res[[i]],syst[[i]][j])
+ res[[i]]<-paste(res[[i]],"+")
+ }
+ res[[i]]<-paste(res[[i]],syst[[i]][length(syst[[i]])])
+ res[[i]]<-as.formula(res[[i]])
+
+ }
+ return (res)
+}
diff --git a/R/systemfit/param.multiple.R b/R/systemfit/param.multiple.R
new file mode 100644
index 0000000..6b26b0a
--- /dev/null
+++ b/R/systemfit/param.multiple.R
@@ -0,0 +1,19 @@
+param.multiple <- function(object, num = NULL, bootstrap = FALSE) {
+ if (!bootstrap) {
+ coef <- mvrnorm(num, mu=coef(object), Sigma=vcov(object))
+ if (object$zelig %in% c("sur","2sls","w2sls","3sls")) {
+ res <- coef
+ }
+ else
+ res <- coef
+ }
+ else {
+ coef <- coef(object)
+ res <- coef
+ }
+ res
+}
+
+
+
+
diff --git a/R/systemfit/plot.zelig.2sls.R b/R/systemfit/plot.zelig.2sls.R
new file mode 100644
index 0000000..973ba9a
--- /dev/null
+++ b/R/systemfit/plot.zelig.2sls.R
@@ -0,0 +1,14 @@
+plot.zelig.2sls <- function(x, xlab = "", user.par = FALSE, ...) {
+ k <- length(x$qi)
+ op <- par(no.readonly = TRUE)
+ if (!user.par)
+ par(mar = c(4,4,2,1), tcl = -0.25, mgp = c(2, 0.6, 0))
+ par(mfrow = c(k, dims(x$qi[[1]])[2]))
+ for (i in 1:k) {
+ for (j in 1:dims(x$qi[[i]])[2]){
+ qi <- as.vector((x$qi[[i]])[,j])
+ plot(density(qi), main = x$qi.name[[i]], xlab = xlab, ...)
+ }
+ }
+ par(op)
+}
diff --git a/R/systemfit/plot.zelig.3sls.R b/R/systemfit/plot.zelig.3sls.R
new file mode 100644
index 0000000..8507059
--- /dev/null
+++ b/R/systemfit/plot.zelig.3sls.R
@@ -0,0 +1,14 @@
+plot.zelig.3sls <- function(x, xlab = "", user.par = FALSE, ...) {
+ k <- length(x$qi)
+ op <- par(no.readonly = TRUE)
+ if (!user.par)
+ par(mar = c(4,4,2,1), tcl = -0.25, mgp = c(2, 0.6, 0))
+ par(mfrow = c(k, dims(x$qi[[1]])[2]))
+ for (i in 1:k) {
+ for (j in 1:dims(x$qi[[i]])[2]){
+ qi <- as.vector((x$qi[[i]])[,j])
+ plot(density(qi), main = x$qi.name[[i]], xlab = xlab, ...)
+ }
+ }
+ par(op)
+}
diff --git a/R/systemfit/plot.zelig.sur.R b/R/systemfit/plot.zelig.sur.R
new file mode 100644
index 0000000..691e8e1
--- /dev/null
+++ b/R/systemfit/plot.zelig.sur.R
@@ -0,0 +1,14 @@
+plot.zelig.sur <- function(x, xlab = "", user.par = FALSE, ...) {
+ k <- length(x$qi)
+ op <- par(no.readonly = TRUE)
+ if (!user.par)
+ par(mar = c(4,4,2,1), tcl = -0.25, mgp = c(2, 0.6, 0))
+ par(mfrow = c(k, dims(x$qi[[1]])[2]))
+ for (i in 1:k) {
+ for (j in 1:dims(x$qi[[i]])[2]){
+ qi <- as.vector((x$qi[[i]])[,j])
+ plot(density(qi), main = x$qi.name[[i]], xlab = xlab, ...)
+ }
+ }
+ par(op)
+}
diff --git a/R/systemfit/plot.zelig.w2sls.R b/R/systemfit/plot.zelig.w2sls.R
new file mode 100644
index 0000000..925a660
--- /dev/null
+++ b/R/systemfit/plot.zelig.w2sls.R
@@ -0,0 +1,14 @@
+plot.zelig.w2sls <- function(x, xlab = "", user.par = FALSE, ...) {
+ k <- length(x$qi)
+ op <- par(no.readonly = TRUE)
+ if (!user.par)
+ par(mar = c(4,4,2,1), tcl = -0.25, mgp = c(2, 0.6, 0))
+ par(mfrow = c(k, dims(x$qi[[1]])[2]))
+ for (i in 1:k) {
+ for (j in 1:dims(x$qi[[i]])[2]){
+ qi <- as.vector((x$qi[[i]])[,j])
+ plot(density(qi), main = x$qi.name[[i]], xlab = xlab, ...)
+ }
+ }
+ par(op)
+}
diff --git a/R/systemfit/qi.multiple.R b/R/systemfit/qi.multiple.R
new file mode 100644
index 0000000..7f69581
--- /dev/null
+++ b/R/systemfit/qi.multiple.R
@@ -0,0 +1,60 @@
+qi.multiple <- function(object, simpar, x, x1 = NULL, y = NULL) {
+
+ check <- FALSE
+ model <- object$zelig
+ coef<-list()
+ om<-attr(object$terms,"omit")
+ nreq<-nrow(om)
+
+ start<-1
+ for(i in 1:nreq){
+ eqni<-paste("eqn",i,sep="")
+ coef[[i]]<-simpar[,start:(start+length(attr(x,eqni))-1)]
+ start<-start+length(attr(x,eqni))
+ }
+
+ fillmatrix<-function(simpar,x,nreq){
+ r<-list()
+ eta<-list()
+ if(nrow(x)==1)
+ q<-array(NA,c(nrow(simpar),nreq))
+ else
+ q<- array(NA,c(nrow(simpar),nreq,nrow(x)))
+
+ for(i in 1:nreq){
+ eqn=paste("eqn",i,sep="")
+ r[[i]]= as.matrix(x[,attr(x,eqn)])
+ if(ncol(r[[i]])==1){
+ eta[[i]] <- coef[[i]] %*% r[[i]]
+ q[,i] <- eta[[i]]
+ }
+ else
+ {
+ eta[[i]] <- coef[[i]] %*% t(r[[i]])
+ q[,i,] <- eta[[i]]
+ }
+ }
+ return (q)
+
+ }
+ pr<-ev<-fillmatrix(simpar,x,nreq)
+
+ qi <- list(ev = ev,pr=pr)
+ qi.name <- list(ev = "Expected Values: E(Y|X)",pr = "Predicted Values: Y|X")
+
+
+ if (!is.null(x1)){
+ theta1<-fillmatrix(simpar,x1,nreq)
+ ev1 <- theta1
+ qi$fd <- ev1-ev
+ qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
+
+ }
+ list(qi=qi, qi.name=qi.name)
+}
+
+
+
+
+
+
diff --git a/R/systemfit/zelig22sls.R b/R/systemfit/zelig22sls.R
new file mode 100644
index 0000000..a745ebf
--- /dev/null
+++ b/R/systemfit/zelig22sls.R
@@ -0,0 +1,15 @@
+zelig22sls <- function(formula, model, data, M,
+ omit = NULL, ...) {
+ check <- library()
+ if(any(check$results[,"Package"] == "systemfit"))
+ require(systemfit)
+ else
+ stop("Please install systemfit using \n install.packages(\"systemfit\")")
+ mf <- match.call(expand.dots = TRUE)
+ mf[[1]] <- as.name("callsystemfit")
+ tmp <- cmsystemfit(formula, omit)
+ mf$eqns <- tmp
+ mf$method<-"2SLS"
+ mf$model<- mf$M<-NULL
+ as.call(mf)
+}
diff --git a/R/systemfit/zelig23sls.R b/R/systemfit/zelig23sls.R
new file mode 100644
index 0000000..5ae0b53
--- /dev/null
+++ b/R/systemfit/zelig23sls.R
@@ -0,0 +1,15 @@
+zelig23sls <- function(formula, model, data, M,
+ omit = NULL, ...) {
+ check <- library()
+ if(any(check$results[,"Package"] == "systemfit"))
+ require(systemfit)
+ else
+ stop("Please install systemfit using \n install.packages(\"systemfit\")")
+ mf <- match.call(expand.dots = TRUE)
+ mf[[1]] <- as.name("callsystemfit")
+ tmp <- cmsystemfit(formula, omit)
+ mf$eqns <- tmp
+ mf$method<-"3SLS"
+ mf$model<- mf$M<-NULL
+ as.call(mf)
+}
diff --git a/R/systemfit/zelig2sur.R b/R/systemfit/zelig2sur.R
new file mode 100644
index 0000000..de353cc
--- /dev/null
+++ b/R/systemfit/zelig2sur.R
@@ -0,0 +1,15 @@
+zelig2sur <- function(formula, model, data, M,
+ omit = NULL,...) {
+ check <- library()
+ if(any(check$results[,"Package"] == "systemfit"))
+ require(systemfit)
+ else
+ stop("Please install systemfit using \n install.packages(\"systemfit\")")
+ mf <- match.call(expand.dots = TRUE)
+ mf[[1]] <- as.name("callsystemfit")
+ tmp <- cmsystemfit(formula, omit)
+ mf$eqns <- tmp
+ mf$method<-"SUR"
+ mf$model<- mf$M<-NULL
+ as.call(mf)
+}
diff --git a/R/systemfit/zelig2w2sls.R b/R/systemfit/zelig2w2sls.R
new file mode 100644
index 0000000..438f507
--- /dev/null
+++ b/R/systemfit/zelig2w2sls.R
@@ -0,0 +1,15 @@
+zelig2w2sls <- function(formula, model, data, M,
+ omit = NULL, ...) {
+ check <- library()
+ if(any(check$results[,"Package"] == "systemfit"))
+ require(systemfit)
+ else
+ stop("Please install systemfit using \n install.packages(\"systemfit\")")
+ mf <- match.call(expand.dots = TRUE)
+ mf[[1]] <- as.name("callsystemfit")
+ tmp <- cmsystemfit(formula, omit)
+ mf$eqns <- tmp
+ mf$method<-"3SLS"
+ mf$model<- mf$M<-NULL
+ as.call(mf)
+}
diff --git a/R/vcov.relogit.R b/R/vcov.relogit.R
index 5598e34..db6bce1 100644
--- a/R/vcov.relogit.R
+++ b/R/vcov.relogit.R
@@ -1,22 +1,2 @@
-vcov.relogit <- function(object, ...) {
- if (is.null(object$call$bias.correct))
- object$call$bias.correct <- TRUE
- dta <- eval(object$call$data, sys.parent())
- tt <- terms(object)
- dta <- dta[complete.cases(model.frame(tt, dta)),]
- n <- nrow(dta)
- k <- ncol(dta)
- bias.cor <- function(res.obj, bc, n, k){
- if (bc) {
- res.obj$cov.unscaled <- res.obj$cov.unscaled * (n/(n+k))^2
- res.obj$cov.scaled <- res.obj$cov.unscaled * res.obj$dispersion
- res.obj$coef[,2] <- sqrt(diag(res.obj$cov.scaled))
- res.obj$coef[,3] <- res.obj$coef[,1] / res.obj$coef[,2]
- res.obj$coef[,4 ] <- 2*pt(-abs(res.obj$coef[,3]), res.obj$df.residual)
- }
- res.obj
- }
- res <- summary.glm(object, ...)
- res <- bias.cor(res, object$call$bias.correct, n, k)
- res$cov.scaled
-}
+vcov.relogit <- function(object, ...)
+ summary.relogit(object, ...)$cov.scaled
diff --git a/R/zelig.R b/R/zelig.R
index d8b90f6..e1140b1 100644
--- a/R/zelig.R
+++ b/R/zelig.R
@@ -46,7 +46,7 @@ zelig <- function(formula, model, data, by = NULL, ...) {
mf$data <- d
res <- eval(as.call(mf))
if (exists(fn2))
- res <- do.call(fn2, list(res = res, fcall = as.list(mf),
+ res <- do.call(fn2, list(res = res, fcall = mf,
zcall = as.list(zelig.call)))
res$call <- zelig.call
res$data <- res$call$data
diff --git a/R/zelig2MCMC.R b/R/zelig2MCMC.R
new file mode 100644
index 0000000..5732e7e
--- /dev/null
+++ b/R/zelig2MCMC.R
@@ -0,0 +1,348 @@
+zelig2ei.hier <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 50000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 5000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf$model <- mf$M <- NULL
+ temp <- mcmcei(formula=formula, data=data)
+
+ if ((any(temp<0)) || ((any(temp<1) && !any(temp==0) ) && any(temp>1)))
+ stop("data need to be either counts or proportions.\n")
+ if (is.null(mf$N))
+ {
+ if (all(temp>=0)) #N is not needed
+ {
+ mf$r0 <- temp$r0
+ mf$r1 <- temp$r1
+ mf$c0 <- temp$c0
+ mf$c1 <- temp$c1
+ }
+ else stop("Needs total counts for inputs as porportion.\n")
+ }
+ else if (((length(mf$N)!= nrow(data)) && (length(mf$N)!=1)) || (any(mf$N<1)))
+ stop("N needs to have same length as the observations and be postive numbers\n.")
+ else if ((all(temp<1)) && (all(mf$N>1)))
+ {
+ mf$r0 <- round(temp$r0*mf$N)
+ mf$r1 <- mf$N-mf$r0
+ mf$c0 <- round(temp$c0*mf$N)
+ mf$c1 <- mf$N-mf$c0
+
+ }
+
+ mf[[1]] <- MCMCpack::MCMChierEI
+ as.call(mf)
+}
+
+zelig2ei.dynamic <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 50000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 5000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf$model <- mf$M <- NULL
+ temp <- mcmcei(formula=formula, data=data)
+
+ if ((any(temp<0)) || ((any(temp<1) && !any(temp==0) ) && any(temp>1)))
+ stop("data need to be either counts or proportions.\n")
+ if (is.null(mf$N))
+ {
+ if (all(temp>=0)) #N is not needed
+ {
+ mf$r0 <- temp$r0
+ mf$r1 <- temp$r1
+ mf$c0 <- temp$c0
+ mf$c1 <- temp$c1
+ }
+ else stop("Needs total counts for inputs as porportion.\n")
+ }
+ else if (((length(mf$N)!= nrow(data)) && (length(mf$N)!=1)) || (any(mf$N<1)))
+ stop("N needs to have same length as the observations and be postive numbers\n.")
+ else if ((all(temp<1)) && (all(mf$N>1)))
+ {
+ mf$r0 <- round(temp$r0*mf$N)
+ mf$r1 <- mf$N-mf$r0
+ mf$c0 <- round(temp$c0*mf$N)
+ mf$c1 <- mf$N-mf$c0
+
+ }
+
+ mf[[1]] <- MCMCpack::MCMCdynamicEI
+ as.call(mf)
+}
+
+
+
+
+zelig2logit.bayes <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf$model <- mf$M <- NULL
+
+ mf[[1]] <- MCMCpack::MCMClogit
+ as.call(mf)
+}
+
+zelig2probit.bayes <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf$model <- mf$M <- NULL
+
+ mf[[1]] <- MCMCpack::MCMCprobit
+ as.call(mf)
+}
+
+zelig2normal.bayes <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+ mf$model <- mf$M <- NULL
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf[[1]] <- MCMCpack::MCMCregress
+ as.call(mf)
+}
+
+zelig2poisson.bayes <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+ mf$model <- mf$M <- NULL
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+ mf[[1]] <- MCMCpack::MCMCpoisson
+ as.call(mf)
+}
+
+zelig2tobit.bayes <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf$model <- mf$M <- NULL
+ mf[[1]] <- MCMCpack::MCMCtobit
+ as.call(mf)
+}
+
+zelig2mlogit.bayes <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ require(stats)
+ mf <- match.call(expand.dots = TRUE)
+ mf$model <- mf$M <- NULL
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf[[1]] <- MCMCpack::MCMCmnl
+ as.call(mf)
+}
+
+zelig2oprobit.bayes <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ require(stats)
+ mf <- match.call(expand.dots = TRUE)
+ mf$model <- mf$M <- NULL
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf[[1]] <- MCMCpack::MCMCoprobit
+ as.call(mf)
+}
+
+
+
+
+
+
+
+zelig2factor.bayes <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 20000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ if (is.null(mf$factors)) mf$factors<-2
+ else if (mf$factors<2) stop("number of factors needs to be at
+ least 2")
+ mf$model <- mf$M <- NULL
+ mf$x <- as.matrix(model.response(model.frame(formula, data=data, na.action=NULL)))
+ mf[[1]] <- MCMCpack::MCMCfactanal
+ as.call(mf)
+}
+
+zelig2factor.ord <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 20000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+ if (is.null(mf$factors)) mf$factors<-2
+ else if (mf$factors<1) stop("number of factors needs to be at
+ least 1")
+ mf$model <- mf$M <- NULL
+ mf$x <- as.matrix(model.response(model.frame(formula, data=data, na.action=NULL)))
+
+ mf[[1]] <- MCMCpack::MCMCordfactanal
+ as.call(mf)
+}
+
+zelig2factor.mix <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ if (is.null(mf$factors)) mf$factors<-2
+ else if (mf$factors<1) stop("number of factors needs to be at
+ least 1")
+ mf$model <- mf$M <- NULL
+
+ var <- model.response(model.frame(formula, data=data,
+ na.action=NULL))
+ varnames <- colnames(var)
+ mf$x <- as.formula(paste("~", paste(varnames, collapse="+")))
+ mf[[1]] <- MCMCpack::MCMCmixfactanal
+
+ as.call(mf)
+}
+
+
+zelig2irt1d <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 20000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+
+ mf$model <- mf$M <- NULL
+ mf$datamatrix <- as.matrix(model.response(model.frame(formula, data=data,
+ na.action=NULL)))
+ mf$datamatrix <- t(mf$datamatrix)
+ mf[[1]] <- MCMCpack::MCMCirt1d
+
+ as.call(mf)
+}
+
+
+zelig2irtkd <- function(formula, model, data, M, ...) {
+ require(MCMCpack)
+ mf <- match.call(expand.dots = TRUE)
+ if (is.null(mf$verbose) || !mf$verbose) mf$verbose <- 0
+ else
+ {
+ if (is.null(mf$mcmc)) mcmc <- 10000
+ else mcmc <- mf$mcmc
+ if (is.null(mf$burnin)) burnin <- 1000
+ else burnin <- mf$burnin
+ mf$verbose <- round((mcmc+burnin)/10)
+ }
+ if (is.null(mf$dimensions)) mf$dimensions <- 1
+
+ mf$model <- mf$M <- NULL
+ mf$datamatrix <- as.matrix(model.response(model.frame(formula, data=data,
+ na.action=NULL)))
+ mf$datamatrix <- t(mf$datamatrix)
+ mf[[1]] <- MCMCpack::MCMCirtKd
+
+ as.call(mf)
+}
+
diff --git a/R/zelig2exp.R b/R/zelig2exp.R
index d8becda..7ff4dca 100644
--- a/R/zelig2exp.R
+++ b/R/zelig2exp.R
@@ -4,12 +4,14 @@ zelig2exp <- function(formula, model, data, M, ...) {
require(survival)
mf[[1]] <- survival::survreg
mf$dist <- "exponential"
+ if (is.null(mf$robust))
+ mf$robust <- FALSE
if (!is.null(mf$cluster) & !mf$robust)
stop("\nIf cluster is specified, robust must be TRUE.")
if (!is.null(mf$cluster)) {
mf$formula <- as.formula(paste(paste(deparse(formula[[2]])),
paste(deparse(formula[[1]])),
- paste(deparse(formula[[3]])),
+ paste(deparse(formula[[3]], width.cutoff=500)),
paste("+", " cluster(",
mf$cluster, ")")))
mf$cluster <- NULL
@@ -17,9 +19,8 @@ zelig2exp <- function(formula, model, data, M, ...) {
else if (mf$robust)
mf$formula <- as.formula(paste(paste(deparse(formula[[2]])),
paste(deparse(formula[[1]])),
- paste(deparse(formula[[3]])),
+ paste(deparse(formula[[3]], width.cutoff=500)),
paste("+", " cluster(1:nrow(",
- deparse(formula[[2]]),
- "))")))
+ deparse(formula[[2]]),"))")))
as.call(mf)
}
diff --git a/R/zelig2logit.R b/R/zelig2logit.R
index d176d79..8585cb8 100644
--- a/R/zelig2logit.R
+++ b/R/zelig2logit.R
@@ -1,7 +1,9 @@
zelig2logit <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
mf$model <- mf$M <- mf$robust <- NULL
- mf[[1]] <- glm
+ mf$formula <- as.formula(paste("cbind(", formula[[2]], ", 1 -", formula[[2]], ")",
+ "~", deparse(formula[[3]], width.cutoff=500), sep = ""))
+ mf[[1]] <- stats::glm
mf$family <- binomial(link="logit")
as.call(mf)
}
diff --git a/R/zelig2lognorm.R b/R/zelig2lognorm.R
index 2c56a0f..9ffd08f 100644
--- a/R/zelig2lognorm.R
+++ b/R/zelig2lognorm.R
@@ -4,12 +4,14 @@ zelig2lognorm <- function(formula, model, data, M, ...) {
require(survival)
mf[[1]] <- survival::survreg
mf$dist <- "lognormal"
+ if (is.null(mf$robust))
+ mf$robust <- FALSE
if (!is.null(mf$cluster) & !mf$robust)
stop("\nIf cluster is specified, robust must be TRUE.")
if (!is.null(mf$cluster)) {
mf$formula <- as.formula(paste(paste(deparse(formula[[2]])),
paste(deparse(formula[[1]])),
- paste(deparse(formula[[3]])),
+ paste(deparse(formula[[3]], width.cutoff=500)),
paste("+", " cluster(",
mf$cluster, ")")))
mf$cluster <- NULL
@@ -17,7 +19,7 @@ zelig2lognorm <- function(formula, model, data, M, ...) {
else if (mf$robust)
mf$formula <- as.formula(paste(paste(deparse(formula[[2]])),
paste(deparse(formula[[1]])),
- paste(deparse(formula[[3]])),
+ paste(deparse(formula[[3]], width.cutoff=500)),
paste("+", " cluster(1:nrow(",
deparse(formula[[2]]),
"))")))
diff --git a/R/zelig2negbin.R b/R/zelig2negbin.R
index c4cff1d..7c64930 100644
--- a/R/zelig2negbin.R
+++ b/R/zelig2negbin.R
@@ -1,6 +1,6 @@
zelig2negbin <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- NULL
+ mf$model <- mf$M <- mf$robust <- NULL
mf[[1]] <- MASS::glm.nb
as.call(mf)
}
diff --git a/R/zelig2probit.R b/R/zelig2probit.R
index 9c35621..e9ddc53 100644
--- a/R/zelig2probit.R
+++ b/R/zelig2probit.R
@@ -1,6 +1,8 @@
zelig2probit <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
mf$model <- mf$M <- mf$probit <- NULL
+ mf$formula <- as.formula(paste("cbind(", formula[[2]], ", 1 -", formula[[2]], ")",
+ "~", deparse(formula[[3]], width.cutoff=500), sep = ""))
mf[[1]] <- stats::glm
mf$family <- binomial(link="probit")
as.call(mf)
diff --git a/R/zelig2relogit.R b/R/zelig2relogit.R
index 93fedc3..116227a 100644
--- a/R/zelig2relogit.R
+++ b/R/zelig2relogit.R
@@ -1,6 +1,10 @@
zelig2relogit <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- NULL
+ mf$model <- mf$M <- mf$robust <- NULL
+ if (is.null(mf$case.control))
+ mf$case.control <- "prior"
+ mf$formula <- as.formula(paste("cbind(", formula[[2]], ", 1 -", formula[[2]], ")",
+ "~", deparse(formula[[3]], width.cutoff=500), sep = ""))
if (is.null(mf$bias.correct))
mf$bias.correct <- TRUE
mf[[1]] <- as.name("relogit")
diff --git a/R/zelig2weibull.R b/R/zelig2weibull.R
index 75e46e1..6b0e23e 100644
--- a/R/zelig2weibull.R
+++ b/R/zelig2weibull.R
@@ -4,12 +4,14 @@ zelig2weibull <- function(formula, model, data, M, ...) {
require(survival)
mf[[1]] <- survival::survreg
mf$dist <- "weibull"
+ if (is.null(mf$robust))
+ mf$robust <- FALSE
if (!is.null(mf$cluster) & !mf$robust)
stop("\nIf cluster is specified, robust must be TRUE.")
if (!is.null(mf$cluster)) {
mf$formula <- as.formula(paste(paste(deparse(formula[[2]])),
paste(deparse(formula[[1]])),
- paste(deparse(formula[[3]])),
+ paste(deparse(formula[[3]], width.cutoff=500)),
paste("+", " cluster(",
mf$cluster, ")")))
mf$cluster <- NULL
@@ -17,7 +19,7 @@ zelig2weibull <- function(formula, model, data, M, ...) {
else if (mf$robust)
mf$formula <- as.formula(paste(paste(deparse(formula[[2]])),
paste(deparse(formula[[1]])),
- paste(deparse(formula[[3]])),
+ paste(deparse(formula[[3]], width.cutoff=500)),
paste("+", " cluster(1:nrow(",
deparse(formula[[2]]),
"))")))
diff --git a/R/zelig3MCMC.R b/R/zelig3MCMC.R
new file mode 100644
index 0000000..5d48c24
--- /dev/null
+++ b/R/zelig3MCMC.R
@@ -0,0 +1,136 @@
+zelig3ei.hier <- zelig3ei.dynamic <- function(res, fcall=NULL, zcall=NULL) {
+
+ out <- list()
+ out$coefficients <- res
+ out$formula <- zcall$formula
+ out$data <- zcall$data
+
+ if (!is.null(zcall$N))
+ out$N <- zcall$N
+
+ out$model <- model.frame(formula=eval(out$formula),
+ data=eval(out$data))
+ out$terms <- attr(out$model, "terms")
+ attr(out$terms,"intercept") <- 0
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
+ class(out) <- "MCMCZelig"
+
+ out
+}
+
+
+zelig3logit.bayes <- zelig3oprobit.bayes <- zelig3poisson.bayes <-
+ zelig3mlogit.bayes <- zelig3normal.bayes <- function(res, fcall=NULL, zcall=NULL) {
+
+ out <- list()
+ out$coefficients <- res
+ out$formula <- zcall$formula
+ out$data <- zcall$data
+
+ out$model <- model.frame(formula=eval(out$formula),
+ data=eval(out$data))
+ out$terms <- attr(out$model, "terms")
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
+
+ class(out) <- "MCMCZelig"
+
+ out
+}
+
+zelig3probit.bayes <- function(res, fcall=NULL, zcall=NULL) {
+
+ out <- list()
+ if (is.null(zcall$bayes.resid))
+ zcall$bayes.resid <- FALSE
+
+ if (zcall$bayes.resid==FALSE)
+ out$coefficients <- res
+ else
+ {
+ p<-dim(model.matrix(eval(zcall$formula), eval(zcall$data)))[2]
+ out$coefficients <- res[,1:p]
+ out$bayes.residuals <- res[, -(1:p)]
+ }
+
+ out$formula <- zcall$formula
+ out$data <- zcall$data
+
+ out$model <- model.frame(formula=eval(out$formula),data=eval(out$data))
+ out$terms <- attr(out$model, "terms")
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
+
+ class(out) <- "MCMCZelig"
+
+ out
+
+ }
+
+ zelig3tobit.bayes <- function(res, fcall=NULL, zcall=NULL) {
+
+ out <- list()
+ out$coefficients <- res
+ out$formula <- zcall$formula
+ if (!is.null(zcall$below)) out$below <- zcall$below
+ else out$below <- 0
+
+ if (!is.null(zcall$above)) out$above <- zcall$above
+ else out$above <- Inf
+
+ out$data <- zcall$data
+
+ out$model <- model.frame(formula=eval(out$formula),
+ data=eval(out$data))
+ out$terms <- attr(out$model, "terms")
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
+
+ class(out) <- "MCMCZelig"
+
+ out
+}
+
+ zelig3factor.bayes <- zelig3factor.ord <- zelig3factor.mix <- function(res, fcall=NULL, zcall=NULL) {
+
+ out <- list()
+ out$coefficients <- res
+ out$formula <- zcall$formula
+ out$data <- zcall$data
+
+ out$model <- model.frame(formula=eval(out$formula),
+ data=eval(out$data))
+ out$terms <- attr(out$model, "terms")
+ attr(out$terms,"intercept") <- 0
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
+
+ class(out) <- "MCMCZelig"
+
+ out
+}
+
+
+ zelig3irt1d <- zelig3irtkd <- function(res, fcall=NULL, zcall=NULL) {
+
+ out <- list()
+ out$coefficients <- res
+ out$formula <- zcall$formula
+ out$data <- zcall$data
+
+ out$model <- model.frame(formula=eval(out$formula),
+ data=eval(out$data))
+ out$terms <- attr(out$model, "terms")
+ attr(out$terms,"intercept") <- 0
+ if (is.null(zcall$seed)) out$seed <- NA
+ else out$seed <- zcall$seed
+
+ class(out) <- "MCMCZelig"
+
+ out
+}
+
+
+
+
diff --git a/R/zelig3relogit.R b/R/zelig3relogit.R
new file mode 100644
index 0000000..3d9a482
--- /dev/null
+++ b/R/zelig3relogit.R
@@ -0,0 +1,53 @@
+zelig3relogit <- function(res, fcall = NULL, zcall = NULL) {
+
+ if ("relogit2" %in% class(res)) {
+ zcall$tau <- res$tau1$tau
+ res$tau1$call <- as.call(zcall)
+ zcall$tau <- res$tau2$tau
+ res$tau2$call <- as.call(zcall)
+ }
+
+ if (is.null(zcall$robust)) {
+ if (is.null(zcall$weighting))
+ rob <- FALSE
+ else if (zcall$weighting) {
+ warning("robust is set to TRUE because weighting is used")
+ rob <- TRUE
+ }
+ else
+ rob <- FALSE
+ }
+ else if (is.logical(zcall$robust))
+ if (!zcall$robust) {
+ rob <- TRUE
+ warning("robust is set to TRUE because weighting is used")
+ }
+ else
+ rob <- zcall$robust
+ if (is.list(rob)) {
+ require(sandwich)
+ if (!any(rob$method %in% c("vcovHAC", "kernHAC", "weave")))
+ stop("such a robust option is not supported")
+ else {
+ if ("relogit2" %in% class(res)) {
+ class(res$tau1) <- class(res$tau2) <- c("relogit", "glm.robust")
+ res$tau1$robust <- res$tau2$robust <- rob
+ }
+ else {
+ class(res) <- c("relogit", "glm.robust")
+ res$robust <- rob
+ }
+ }
+ }
+ else if (!is.logical(rob))
+ stop("invalid input for robust. Choose either TRUE or a list of options.")
+ else if (rob) {
+ require(sandwich)
+ if ("relogit2" %in% class(res))
+ class(res$tau1) <- class(res$tau2) <- c("relogit", "glm.robust")
+ else
+ class(res) <- c("relogit", "glm.robust")
+ }
+ return(res)
+}
+
diff --git a/README b/README
index e974adc..970c5db 100644
--- a/README
+++ b/README
@@ -1,3 +1,25 @@
+2.4-1 (August 15, 2005): Stable release for R 2.0.0-2.1.1. Added the
+ following Bayesian models: factor analysis, mixed factor analysis,
+ ordinal factor analysis, unidimensional item response theory,
+ k-dimensional item response theory, logit, multinomial logit, normal,
+ ordinal probit, Poisson, and tobit. Also fixed minor bug in
+ formula (long variable names coerced to list).
+
+2.3-2 (August 5, 2005): Stable release for R 2.0.0-2.1.1. Fixed bug in
+ simulation procedure for lognormal model.
+
+2.3-1 (August 4, 2005): Stable release for R 2.0.0-2.1.1.
+ Fixed documentation errors related to model parameterization and code
+ bugs related to first differences and conditional prediction for
+ exponential, lognormal, and Weibull models. (reported by Alison Post)
+
+2.2-4 (July 30, 2005): Stable release for R 2.0.0-2.1.1. Revised
+ relogit, additing options for weightig in addition to prior
+ correction. (reported by Martin Plöderl)
+
+2.2-3 (July 24, 2005): Stable release for R 2.0.0-2.1.1. Fixed bug
+ associated with robust standard errors for negative binomial.
+
2.2-2 (July 13, 2005): Stable release for R 2.0.0-2.1.1. Fixed bug in
setx(). (reported by Ying Lu)
diff --git a/data/PErisk.txt b/data/PErisk.txt
new file mode 100644
index 0000000..2e556ca
--- /dev/null
+++ b/data/PErisk.txt
@@ -0,0 +1,63 @@
+"country" "courts" "barb2" "prsexp2" "prscorr2" "gdpw2"
+"Argentina" "Argentina" "0" -0.7207754 "1" "3" 9.69017
+"Australia" "Australia" "1" -6.907755 "5" "4" 10.30484
+"Austria" "Austria" "1" -4.910337 "5" "4" 10.10094
+"Bangladesh" "Bangladesh" "0" 0.7759748 "1" "0" 8.379768
+"Belgium" "Belgium" "1" -4.617344 "5" "4" 10.25012
+"Bolivia" "Bolivia" "0" -2.46144 "0" "0" 8.583543
+"Botswana" "Botswana" "1" -1.244868 "4" "3" 8.77771
+"Brazil" "Brazil" "1" -0.4570337 "4" "3" 9.375601
+"Burma" "Burma" "0" 1.604343 "3" "1" 7.096721
+"Cameroon" "Cameroon" "0" -4.229065 "3" "1" 8.120886
+"Canada" "Canada" "1" -6.907755 "5" "5" 10.41018
+"Chile" "Chile" "1" -1.542761 "3" "2" 9.261224
+"Colombia" "Colombia" "0" -2.057821 "3" "2" 9.191973
+"Congo-Kinshasa" "Congo-Kinshasa" "0" -2.323288 "1" "0" 7.095064
+"Costa Rica" "Costa Rica" "1" -5.090003 "3" "4" 9.167329
+"Cote d'Ivoire" "Cote d'Ivoire" "1" -4.229065 "4" "2" 8.228711
+"Denmark" "Denmark" "1" -6.907755 "5" "5" 10.10651
+"Dominican Republic" "Dominican Republic" "0" -2.378862 "2" "2" 8.899731
+"Ecuador" "Ecuador" "1" -1.845337 "3" "2" 9.117786
+"Finland" "Finland" "1" -6.907755 "5" "5" 10.12367
+"Gambia, The" "Gambia, The" "0" -1.543332 "4" "2" 7.501082
+"Ghana" "Ghana" "0" -1.011517 "2" "1" 7.597396
+"Greece" "Greece" "0" -2.073237 "3" "3" 9.701494
+"Hungary" "Hungary" "1" -0.9041942 "4" "3" 9.35184
+"India" "India" "0" -2.105104 "4" "2" 7.970049
+"Indonesia" "Indonesia" "0" -2.100232 "3" "0" 8.39231
+"Iran" "Iran" "0" 2.337425 "0" "2" 9.368114
+"Ireland" "Ireland" "1" -6.907755 "5" "4" 9.891465
+"Israel" "Israel" "0" -2.319996 "4" "4" 10.06777
+"Italy" "Italy" "1" -6.907755 "4" "3" 10.26078
+"Japan" "Japan" "1" -6.907755 "5" "4" 9.892022
+"Kenya" "Kenya" "0" -2.327605 "2" "2" 7.619724
+"Korea, South" "Korea, South" "0" -2.655795 "4" "1" 9.422787
+"Malawi" "Malawi" "0" -1.469424 "3" "3" 7.029973
+"Malaysia" "Malaysia" "1" -3.927949 "4" "3" 9.178953
+"Mexico" "Mexico" "0" -1.657935 "2" "2" 9.661735
+"Morocco" "Morocco" "0" -3.156958 "3" "1" 8.78048
+"New Zealand" "New Zealand" "1" -6.907755 "5" "5" 10.17626
+"Nigeria" "Nigeria" "0" 0.3001068 "1" "1" 7.68708
+"Norway" "Norway" "1" -6.907755 "5" "5" 10.29833
+"Papua New Guinea" "Papua New Guinea" "1" -2.636158 "4" "2" 8.126518
+"Paraguay" "Paraguay" "0" -0.9707628 "3" "0" 8.727616
+"Philippines" "Philippines" "0" -2.964776 "1" "1" 8.384804
+"Poland" "Poland" "1" 1.317021 "3" "3" 9.0524
+"Portugal" "Portugal" "1" -2.459625 "4" "3" 9.444543
+"Sierra Leone" "Sierra Leone" "0" 1.406406 "3" "1" 7.759614
+"Singapore" "Singapore" "1" -4.848516 "5" "5" 9.882724
+"South Africa" "South Africa" "0" -2.17582 "3" "4" 9.191871
+"Spain" "Spain" "1" -6.907755 "5" "3" 10.04733
+"Sri Lanka" "Sri Lanka" "0" -1.864343 "2" "2" 8.627661
+"Sweden" "Sweden" "1" -6.907755 "4" "5" 10.22434
+"Switzerland" "Switzerland" "1" -6.907755 "5" "5" 10.3411
+"Syria" "Syria" "0" 1.725166 "1" "1" 9.664151
+"Thailand" "Thailand" "0" -6.907755 "3" "2" 8.548692
+"Togo" "Togo" "0" -4.229065 "4" "1" 7.331715
+"Tunisia" "Tunisia" "0" -2.585399 "2" "2" 9.047586
+"Turkey" "Turkey" "0" -2.673243 "3" "2" 8.978912
+"United Kingdom" "United Kingdom" "1" -6.907755 "5" "5" 10.12727
+"Uruguay" "Uruguay" "0" -2.127775 "2" "2" 9.414342
+"Venezuela" "Venezuela" "1" 0.428845 "3" "2" 9.84882
+"Zambia" "Zambia" "0" 0.9658105 "3" "1" 7.726213
+"Zimbabwe" "Zimbabwe" "0" -0.6403214 "3" "2" 7.965893
diff --git a/data/SupremeCourt.txt b/data/SupremeCourt.txt
new file mode 100644
index 0000000..07b11a7
--- /dev/null
+++ b/data/SupremeCourt.txt
@@ -0,0 +1,44 @@
+"Rehnquist" "Stevens" "OConnor" "Scalia" "Kennedy" "Souter" "Thomas" "Ginsburg" "Breyer"
+"1" 0 1 1 0 1 1 0 1 1
+"2" 0 1 0 0 0 1 0 1 1
+"3" 0 1 0 0 0 1 0 1 1
+"4" 0 0 0 0 0 0 0 0 1
+"5" 1 1 0 0 1 0 0 0 0
+"6" 0 1 0 0 0 0 0 0 0
+"7" 0 1 1 0 0 1 0 1 1
+"8" 0 1 0 0 0 0 0 0 0
+"9" 0 1 0 0 1 1 0 1 1
+"10" 1 1 1 0 1 1 0 1 1
+"11" 0 1 1 0 1 1 0 1 1
+"12" 0 1 0 0 0 1 0 1 1
+"13" 1 0 1 1 1 1 1 1 0
+"14" 0 1 0 0 0 1 0 1 1
+"15" 0 1 1 0 0 1 0 1 1
+"16" 0 1 0 0 0 1 0 1 1
+"17" 0 1 1 0 0 0 0 1 1
+"18" 0 1 0 0 0 1 0 1 1
+"19" 0 1 0 0 0 1 0 1 1
+"20" 0 0 0 0 0 0 0 1 0
+"21" 0 1 0 1 0 0 1 0 1
+"22" 0 1 1 0 1 1 0 1 1
+"23" 1 NA NA 0 1 1 0 1 1
+"24" 0 1 0 0 0 1 0 1 1
+"25" 1 1 1 0 1 1 0 1 1
+"26" 0 1 0 0 0 1 0 1 1
+"27" 0 1 1 0 1 1 0 1 1
+"28" 0 1 0 0 0 0 0 0 0
+"29" 0 0 0 1 0 1 1 1 1
+"30" 0 0 1 0 0 1 0 1 1
+"31" 0 1 0 0 0 1 0 1 0
+"32" 0 0 0 0 0 0 0 1 1
+"33" 1 1 1 0 1 1 1 1 1
+"34" 0 1 1 0 0 1 0 1 1
+"35" 0 1 0 0 1 1 0 1 1
+"36" 0 1 0 0 1 1 0 1 1
+"37" 1 1 0 1 1 1 1 0 0
+"38" 0 1 1 0 0 1 0 1 1
+"39" 1 0 1 1 1 1 1 1 0
+"40" 1 0 1 1 1 0 1 0 0
+"41" 0 1 0 0 0 1 0 1 1
+"42" 0 1 0 0 0 1 0 1 1
+"43" 0 1 1 0 0 1 0 1 1
diff --git a/data/eidat.txt b/data/eidat.txt
new file mode 100644
index 0000000..dba9426
--- /dev/null
+++ b/data/eidat.txt
@@ -0,0 +1,11 @@
+"x0" "x1" "t0" "t1"
+"1" 200 3911 2850 1261
+"2" 162 2636 1541 1257
+"3" 206 2460 1091 1575
+"4" 213 1654 517 1350
+"5" 209 637 163 683
+"6" 190 1911 216 1885
+"7" 206 3460 226 3440
+"8" 190 715 102 803
+"9" 183 2058 126 2115
+"10" 189 2658 138 2709
diff --git a/data/newpainters.txt b/data/newpainters.txt
new file mode 100644
index 0000000..373f825
--- /dev/null
+++ b/data/newpainters.txt
@@ -0,0 +1,55 @@
+"Composition" "Drawing" "Colour" "Expression"
+"Da Udine" 200 100 400 100
+"Da Vinci" 400 400 100 400
+"Del Piombo" 100 200 400 300
+"Del Sarto" 200 400 200 300
+"Fr. Penni" 100 400 200 100
+"Guilio Romano" 400 400 100 400
+"Michelangelo" 100 400 100 300
+"Perino del Vaga" 400 400 100 300
+"Perugino" 100 200 300 200
+"Raphael" 400 400 300 400
+"F. Zucarro" 200 200 200 300
+"Fr. Salviata" 300 400 200 300
+"Parmigiano" 200 400 100 300
+"Primaticcio" 400 300 100 300
+"T. Zucarro" 300 300 300 300
+"Volterra" 200 400 100 300
+"Barocci" 300 400 100 300
+"Cortona" 400 300 300 300
+"Josepin" 200 200 100 100
+"L. Jordaens" 300 200 200 300
+"Testa" 200 400 100 300
+"Vanius" 400 400 300 400
+"Bassano" 100 100 400 100
+"Bellini" 100 100 300 100
+"Giorgione" 100 100 400 200
+"Murillo" 100 100 300 200
+"Palma Giovane" 200 100 300 300
+"Palma Vecchio" 100 100 400 100
+"Pordenone" 100 300 400 200
+"Tintoretto" 400 300 400 200
+"Titian" 200 400 400 300
+"Veronese" 400 200 400 100
+"Albani" 300 300 300 300
+"Caravaggio" 100 100 400 100
+"Corregio" 300 200 300 400
+"Domenichino" 400 400 200 400
+"Guercino" 400 200 300 200
+"Lanfranco" 300 200 300 200
+"The Carraci" 400 400 300 400
+"Durer" 100 200 300 300
+"Holbein" 200 200 400 400
+"Pourbus" 100 400 100 300
+"Van Leyden" 100 100 100 200
+"Diepenbeck" 200 200 300 300
+"J. Jordaens" 200 100 400 300
+"Otho Venius" 300 300 300 300
+"Rembrandt" 400 100 400 400
+"Rubens" 400 200 400 400
+"Teniers" 400 200 300 300
+"Van Dyck" 400 200 400 400
+"Bourdon" 200 100 200 200
+"Le Brun" 400 400 200 400
+"Le Suer" 400 400 100 400
+"Poussin" 400 400 100 400
diff --git a/data/swiss.txt b/data/swiss.txt
new file mode 100644
index 0000000..3e8c981
--- /dev/null
+++ b/data/swiss.txt
@@ -0,0 +1,48 @@
+"Fertility" "Agriculture" "Examination" "Education" "Catholic" "Infant.Mortality"
+"Courtelary" 80.2 17 15 12 9.96 22.2
+"Delemont" 83.1 45.1 6 9 84.84 22.2
+"Franches-Mnt" 92.5 39.7 5 5 93.4 20.2
+"Moutier" 85.8 36.5 12 7 33.77 20.3
+"Neuveville" 76.9 43.5 17 15 5.16 20.6
+"Porrentruy" 76.1 35.3 9 7 90.57 26.6
+"Broye" 83.8 70.2 16 7 92.85 23.6
+"Glane" 92.4 67.8 14 8 97.16 24.9
+"Gruyere" 82.4 53.3 12 7 97.67 21
+"Sarine" 82.9 45.2 16 13 91.38 24.4
+"Veveyse" 87.1 64.5 14 6 98.61 24.5
+"Aigle" 64.1 62 21 12 8.52 16.5
+"Aubonne" 66.9 67.5 14 7 2.27 19.1
+"Avenches" 68.9 60.7 19 12 4.43 22.7
+"Cossonay" 61.7 69.3 22 5 2.82 18.7
+"Echallens" 68.3 72.6 18 2 24.2 21.2
+"Grandson" 71.7 34 17 8 3.3 20
+"Lausanne" 55.7 19.4 26 28 12.11 20.2
+"La Vallee" 54.3 15.2 31 20 2.15 10.8
+"Lavaux" 65.1 73 19 9 2.84 20
+"Morges" 65.5 59.8 22 10 5.23 18
+"Moudon" 65 55.1 14 3 4.52 22.4
+"Nyone" 56.6 50.9 22 12 15.14 16.7
+"Orbe" 57.4 54.1 20 6 4.2 15.3
+"Oron" 72.5 71.2 12 1 2.4 21
+"Payerne" 74.2 58.1 14 8 5.23 23.8
+"Paysd'enhaut" 72 63.5 6 3 2.56 18
+"Rolle" 60.5 60.8 16 10 7.72 16.3
+"Vevey" 58.3 26.8 25 19 18.46 20.9
+"Yverdon" 65.4 49.5 15 8 6.1 22.5
+"Conthey" 75.5 85.9 3 2 99.71 15.1
+"Entremont" 69.3 84.9 7 6 99.68 19.8
+"Herens" 77.3 89.7 5 2 100 18.3
+"Martigwy" 70.5 78.2 12 6 98.96 19.4
+"Monthey" 79.4 64.9 7 3 98.22 20.2
+"St Maurice" 65 75.9 9 9 99.06 17.8
+"Sierre" 92.2 84.6 3 3 99.46 16.3
+"Sion" 79.3 63.1 13 13 96.83 18.1
+"Boudry" 70.4 38.4 26 12 5.62 20.3
+"La Chauxdfnd" 65.7 7.7 29 11 13.79 20.5
+"Le Locle" 72.7 16.7 22 13 11.22 18.9
+"Neuchatel" 64.4 17.6 35 32 16.92 23
+"Val de Ruz" 77.6 37.6 15 7 4.97 20
+"ValdeTravers" 67.6 18.7 25 7 8.65 19.5
+"V. De Geneve" 35 1.2 37 53 42.34 18
+"Rive Droite" 44.7 46.6 16 29 50.43 18.2
+"Rive Gauche" 42.8 27.7 22 29 58.33 19.3
diff --git a/data/tobin.txt b/data/tobin.txt
new file mode 100644
index 0000000..f7ce8db
--- /dev/null
+++ b/data/tobin.txt
@@ -0,0 +1,21 @@
+"durable" "age" "quant"
+"1" 0 57.7 236
+"2" 0.7 50.9 283
+"3" 0 48.5 207
+"4" 0 41.7 220
+"5" 0 47.7 238
+"6" 0 59.8 216
+"7" 0 44.3 284
+"8" 3.7 45.1 221
+"9" 0 51.7 275
+"10" 3 50 269
+"11" 10.4 46.8 207
+"12" 0 58 249
+"13" 0 58.9 246
+"14" 0 40 277
+"15" 1.5 34.1 231
+"16" 0 39.9 219
+"17" 0 33.4 240
+"18" 3.5 48.1 266
+"19" 6.1 46.6 214
+"20" 0 53.1 251
diff --git a/demo/00Index b/demo/00Index
index 280023d..f172183 100644
--- a/demo/00Index
+++ b/demo/00Index
@@ -1,24 +1,44 @@
-beta Beta regression and simulation
-blogit Bivariate Logit regression and simulation
-bprobit Bivariate Probit regression and simulation
-conditional Conditional prediction
-exp Exponential regression and simulation
-gamma Gamma regression and simulation
-logit Logit regression and simulation
-lognorm Lognormal regression and simulation
-ls Least Squares regression and simulation
-match Regression and simulation on a matched data set
-mi Regression and simulation on multiply imputed data sets
-mlogit Multinomial Logit regression and simulation
-negbin Negative Binomial regression and simulation
-normal Normal (Gaussian) regression and simulation
-ologit Ordinal Logit regression and simulation
-oprobit Ordinal Probit regression and simulation
-poisson Poisson regression and simulation
-probit Probit regression and simulation
-relogit Rare events logit regression and simulation
-robust Robust estimation and simulation
-strata Regression and simulation in a stratified model
-roc ROC plot
-vertci Vertical confidence interval plot
-weibull Weibull regression and simulation
+beta Beta regression and simulation
+blogit Bivariate Logit regression and simulation
+bprobit Bivariate Probit regression and simulation
+conditional Conditional prediction
+exp Exponential regression and simulation
+gamma Gamma regression and simulation
+logit Logit regression and simulation
+lognorm Lognormal regression and simulation
+ls Least Squares regression and simulation
+match Regression and simulation on a matched data set
+mi Regression and simulation on multiply imputed data sets
+mlogit Multinomial Logit regression and simulation
+negbin Negative Binomial regression and simulation
+normal Normal (Gaussian) regression and simulation
+ologit Ordinal Logit regression and simulation
+oprobit Ordinal Probit regression and simulation
+poisson Poisson regression and simulation
+probit Probit regression and simulation
+relogit Rare events logit regression and simulation
+robust Robust estimation and simulation
+strata Regression and simulation in a stratified model
+roc ROC plot
+vertci Vertical confidence interval plot
+weibull Weibull regression and simulation
+ei.hier Hierarchical Ecological Inference model and simulation
+ei.dynamic Dynamic Ecological Inference model and simulation
+normal.bayes MCMC regression model and simulation
+logit.bayes MCMC logistic regression model and simulation
+probit.bayes MCMC probit regression model and simulation
+tobit.bayes MCMC tobit regression model and simulation
+poisson.bayes MCMC poisson regression model and simulation
+mlogit.bayes MCMC multinomial regression model and simulation
+oprobit.bayes MCMC ordered probit regression model and simulation
+factor.bayes MCMC factor analysis
+factor.ord MCMC factor analysis for ordinal data
+factor.mix MCMC factor analysis for mixed data
+irt1d MCMC one-dimensional item response theory model
+irtkd MCMC K-dimensional item response theory model
+
+
+
+
+
+
diff --git a/demo/ei.dynamic.R b/demo/ei.dynamic.R
new file mode 100644
index 0000000..23334eb
--- /dev/null
+++ b/demo/ei.dynamic.R
@@ -0,0 +1,55 @@
+## Attaching the example dataset:
+data(eidat)
+
+## Estimating the model using MCMCdynamicEI:
+z.out <- zelig(cbind(t0, t1) ~ x0 + x1, model = "ei.dynamic", data = eidat,
+ mcmc = 40000, thin = 10, burnin = 10000, verbose = TRUE)
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+heidel.diag(z.out$coefficients)
+user.prompt()
+raftery.diag(z.out$coefficients)
+user.prompt()
+
+## summarizing the output
+summary(z.out)
+user.prompt()
+
+## Setting values for in-sample simulations given
+## marginal values of X0, X1, T0 and T1:
+x.out <- setx(z.out, fn = NULL, cond=TRUE)
+user.prompt()
+
+## In-sample simulations from the posterior distribution:
+s.out <- sim(z.out, x=x.out)
+
+## Summarizing in-sample simulations at aggregate level
+## weighted by the count in each unit:
+summary(s.out)
+user.prompt()
+## Summarizing in-sample simulations at unit level
+## for the first 5 units:
+summary(s.out, subset = 1:5)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/demo/ei.hier.R b/demo/ei.hier.R
new file mode 100644
index 0000000..c367fce
--- /dev/null
+++ b/demo/ei.hier.R
@@ -0,0 +1,56 @@
+## Attaching the example dataset:
+data(eidat)
+
+## Estimating the model using MCMChierEI:
+z.out <- zelig(cbind(t0, t1) ~ x0 + x1, model="ei.hier", data = eidat,
+ mcmc = 40000, thin = 10, burnin = 10000, verbose = TRUE)
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+heidel.diag(z.out$coefficients)
+user.prompt()
+raftery.diag(z.out$coefficients)
+user.prompt()
+
+## summarizing the output
+summary(z.out)
+user.prompt()
+
+## Setting values for in-sample simulations given
+## marginal values of X0, X1, T0 and T1:
+x.out <- setx(z.out, fn = NULL, cond = TRUE)
+user.prompt()
+
+## In-sample simulations from the posterior distribution:
+s.out <- sim(z.out, x = x.out)
+
+## Summarizing in-sample simulations at aggregate level
+## weighted by the count in each unit:
+summary(s.out)
+user.prompt()
+
+## Summarizing in-sample simulations at unit level
+## for the first 5 units:
+summary(s.out, subset = 1:5)
+user.prompt()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/demo/factor.bayes.R b/demo/factor.bayes.R
new file mode 100644
index 0000000..de85df7
--- /dev/null
+++ b/demo/factor.bayes.R
@@ -0,0 +1,35 @@
+## Attaching the example dataset:
+data(swiss)
+names(swiss) <- c("Fert","Agr","Exam","Educ","Cath","InfMort")
+
+user.prompt()
+
+## Estimating the model using MCMCfactanal:
+z.out <- zelig(cbind(Agr,Exam,Educ,Cath,InfMort)~NULL,
+ model="factor.bayes",
+ data=swiss, factors=2,
+ lambda.constraints=list(Exam=list(1,"+"),
+ Exam=list(2,"-"), Educ=c(2,0),
+ InfMort=c(1,0)),
+ verbose=TRUE, a0=1, b0=0.15,
+ burnin=5000, mcmc=50000)
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+heidel.diag(z.out$coefficients)
+user.prompt()
+raftery.diag(z.out$coefficients)
+user.prompt()
+
+## summarizing the output
+summary(z.out)
+
+
+
+
+
+
+
+
diff --git a/demo/factor.mix.R b/demo/factor.mix.R
new file mode 100644
index 0000000..65964fd
--- /dev/null
+++ b/demo/factor.mix.R
@@ -0,0 +1,26 @@
+## Attaching the example dataset:
+data(PErisk)
+
+## Estimating the model using factor.mix:
+z.out <- zelig(cbind(courts,barb2,prsexp2,prscorr2,gdpw2)~NULL,
+ data=PErisk, model="factor.mix",factors=1,
+ burnin=5000,mcmc=100000, thin=50, verbose=TRUE, L0=0.25)
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+heidel.diag(z.out$coefficients)
+user.prompt()
+
+
+## summarizing the output
+summary(z.out)
+
+
+
+
+
+
+
+
diff --git a/demo/factor.ord.R b/demo/factor.ord.R
new file mode 100644
index 0000000..0bc8321
--- /dev/null
+++ b/demo/factor.ord.R
@@ -0,0 +1,28 @@
+## Attaching the example dataset:
+data(newpainters)
+
+## Estimating the model using factor.ord:
+z.out <- zelig(cbind(Composition,Drawing,Colour,Expression)~NULL,
+ data=newpainters, model="factor.ord",
+ factors=1,
+ burin=1000,mcmc=10000, thin=5, verbose=TRUE,
+ L0=0.5)
+
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+
+heidel.diag(z.out$coefficients)
+user.prompt()
+
+## summarizing the output
+summary(z.out)
+
+
+
+
+
+
+
diff --git a/demo/irt1d.R b/demo/irt1d.R
new file mode 100644
index 0000000..09e898a
--- /dev/null
+++ b/demo/irt1d.R
@@ -0,0 +1,32 @@
+## Attaching the example dataset:
+data(SupremeCourt)
+names(SupremeCourt) <- c("Rehnquist","Stevens","OConnor","Scalia",
+ "Kennedy","Souter","Thomas","Ginsburg","Breyer")
+
+user.prompt()
+
+## Estimating the model using MCMCirt1d:
+z.out <- zelig(cbind(Rehnquist,Stevens,OConnor,Scalia,
+ Kennedy,Souter,Thomas,Ginsburg,Breyer)~NULL,
+ data=SupremeCourt, model="irt1d",
+ B0.alpha=0.2, B0.beta=0.2, burnin=500, mcmc=10000,
+ thin=20, verbose=TRUE)
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+
+heidel.diag(z.out$coefficients)
+user.prompt()
+
+
+## summarizing the output
+summary(z.out)
+
+
+
+
+
+
+
diff --git a/demo/irtkd.R b/demo/irtkd.R
new file mode 100644
index 0000000..d3f6493
--- /dev/null
+++ b/demo/irtkd.R
@@ -0,0 +1,34 @@
+## Attaching the example dataset:
+data(SupremeCourt)
+names(SupremeCourt) <- c("Rehnquist","Stevens","OConnor","Scalia",
+ "Kennedy","Souter","Thomas","Ginsburg","Breyer")
+user.prompt()
+
+## Estimating the model using MCMCirtKd:
+z.out <- zelig(cbind(Rehnquist,Stevens,OConnor,Scalia,
+ Kennedy,Souter,Thomas,Ginsburg,Breyer)~NULL,
+ dimensions=1, data=SupremeCourt, model="irtkd",
+ B0=0.25, burnin=5000, mcmc=50000, thin=10, verbose=TRUE)
+
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+
+heidel.diag(z.out$coefficients)
+user.prompt()
+
+raftery.diag(z.out$coefficients)
+user.prompt()
+
+## summarizing the output
+summary(z.out)
+
+
+
+
+
+
+
+
diff --git a/demo/logit.bayes.R b/demo/logit.bayes.R
new file mode 100644
index 0000000..f1a3a2a
--- /dev/null
+++ b/demo/logit.bayes.R
@@ -0,0 +1,57 @@
+## Attaching the example dataset:
+data(turnout)
+
+## Estimating the model using MCMClogit:
+z.out <- zelig(vote ~ race + educate, model = "logit.bayes",
+ data = turnout, verbose=TRUE)
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+heidel.diag(z.out$coefficients)
+user.prompt()
+raftery.diag(z.out$coefficients)
+user.prompt()
+
+## summarizing the output
+summary(z.out)
+user.prompt()
+
+## Setting values for the explanatory variables to
+## their sample averages:
+x.out <- setx(z.out)
+user.prompt()
+
+## Simulating quantities of interest from the posterior
+## distribution given x.out:
+s.out1 <- sim(z.out, x = x.out)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out1)
+user.prompt()
+
+## Simulating First Differences:
+## Setting education is set to be between low(25th percentile)
+## versus high(75th percentile) while all the other variables
+## held at their default values.
+x.high <- setx(z.out, educate = quantile(turnout$educate, prob = 0.75))
+x.low <- setx(z.out, educate = quantile(turnout$educate, prob = 0.25))
+user.prompt()
+
+## Estimating the first difference for the effect of
+## high versus low trade on unemployment rate:
+s.out2 <- sim(z.out, x = x.high, x1 = x.low)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out2)
+
+
+
+
+
+
+
+
diff --git a/demo/match.R b/demo/match.R
index 788b8db..f6c94f2 100644
--- a/demo/match.R
+++ b/demo/match.R
@@ -2,14 +2,16 @@ library(MatchIt)
data(lalonde)
user.prompt()
## an example for propensity score matching
-match.out1 <- matchit(treat ~ age + educ + black + hispan + married +
- nodegree + re74 + re75, data = lalonde)
+m.out1 <- matchit(treat ~ age + educ + black + hispan + married +
+ nodegree + re74 + re75, data = lalonde)
user.prompt()
-z.out1 <- zelig(re78 ~ pscore, data = match.data(match.out1, "control"), model = "ls")
+z.out1 <- zelig(re78 ~ distance + age + educ + black + hispan + married +
+ nodegree + re74 + re75,
+ data = match.data(m.out1, "control"), model = "ls")
user.prompt()
-x.out1 <- setx(z.out1, fn = NULL, data = match.data(match.out1, "treat"), cond = TRUE)
+x.out1 <- setx(z.out1, fn = NULL, data = match.data(m.out1, "treat"), cond = TRUE)
user.prompt()
s.out1 <- sim(z.out1, x = x.out1)
@@ -18,16 +20,17 @@ summary(s.out1)
user.prompt()
## an example for subclassification
-match.out2 <- matchit(treat ~ age + educ + black + hispan + married +
- nodegree + re74 + re75, data = lalonde,
- replace = FALSE, subclass = 3)
+m.out2 <- matchit(treat ~ age + educ + black + hispan + married +
+ nodegree + re74 + re75, data = lalonde,
+ replace = FALSE, subclass = 3)
user.prompt()
-z.out2 <- zelig(re78 ~ pscore, data = match.data(match.out2, "control"),
- model="ls", by="psclass")
+z.out2 <- zelig(re78 ~ distance + age + educ +
+ married + nodegree + re74 + re75, data = match.data(m.out2, "control"),
+ model="ls", by="subclass")
user.prompt()
-x.out2 <- setx(z.out2, fn = NULL, data = match.data(match.out2, "treat"), cond = TRUE)
+x.out2 <- setx(z.out2, fn = NULL, data = match.data(m.out2, "treat"), cond = TRUE)
user.prompt()
s.out2 <- sim(z.out2, x = x.out2)
diff --git a/demo/mlogit.bayes.R b/demo/mlogit.bayes.R
new file mode 100644
index 0000000..aeb785f
--- /dev/null
+++ b/demo/mlogit.bayes.R
@@ -0,0 +1,58 @@
+## Attaching the example dataset:
+data(mexico)
+
+## Estimating the model using mlogit.bayes:
+z.out <- zelig(vote88 ~ pristr + othcok + othsocok, model = "mlogit.bayes",
+ data = mexico)
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+heidel.diag(z.out$coefficients)
+user.prompt()
+
+raftery.diag(z.out$coefficients)
+user.prompt()
+
+## Summarizing the output
+summary(z.out)
+user.prompt()
+
+## Setting values for the explanatory variables to
+## their sample averages:
+x.out <- setx(z.out)
+user.prompt()
+
+## Simulating quantities of interest from the posterior
+## distribution given x.out:
+s.out1 <- sim(z.out, x = x.out)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out1)
+user.prompt()
+
+## Simulating First Differences:
+## Setting explanatory variables to their default(mean/mode)
+## values, with pristr(the strength of PRI) equal to 1(weak) or
+## 3(strong)
+x.weak <- setx(z.out, pristr = 1)
+x.strong <- setx(z.out, pristr = 3)
+
+user.prompt()
+
+## Estimating the first difference for the effect of
+## military action on the probabilities of
+## incurring differnt level of cost:
+s.out2 <- sim(z.out, x = x.strong, x1 = x.weak)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out2)
+
+
+
+
+
+
+
+
diff --git a/demo/normal.R b/demo/normal.R
index 42de1e0..b10e6d2 100644
--- a/demo/normal.R
+++ b/demo/normal.R
@@ -32,24 +32,24 @@ plot(s.out1)
# Note that you do not need to create dummy variables, as the program
# will automatically parse the unique values in the selected variables
# into dummy variables.
-user.prompt()
-z.out2 <- zelig(unem ~ gdp + trade + capmob + as.factor(year)
- + as.factor(country), model = "normal", data = macro)
+#user.prompt()
+#z.out2 <- zelig(unem ~ gdp + trade + capmob + as.factor(year)
+# + as.factor(country), model = "normal", data = macro)
# Set values for the explanatory variables, using the default mean/mode
# values, with country set to the United States and Japan, respectively:
-user.prompt()
-x.US <- setx(z.out2, country = "United States")
-x.Japan <- setx(z.out2, country = "Japan")
+#user.prompt()
+#x.US <- setx(z.out2, country = "United States")
+#x.Japan <- setx(z.out2, country = "Japan")
# Simulate quantities of interest:
-user.prompt()
-s.out2 <- sim(z.out2, x = x.US, x1 = x.Japan)
-user.prompt()
-summary(s.out2)
+#user.prompt()
+#s.out2 <- sim(z.out2, x = x.US, x1 = x.Japan)
+#user.prompt()
+#summary(s.out2)
# Plot differences:
-user.prompt()
-plot(s.out2)
+#user.prompt()
+#plot(s.out2)
diff --git a/demo/normal.bayes.R b/demo/normal.bayes.R
new file mode 100644
index 0000000..7745f15
--- /dev/null
+++ b/demo/normal.bayes.R
@@ -0,0 +1,58 @@
+## Attaching the example dataset:
+data(macro)
+
+## Estimating the model using normal.bayes:
+z.out <- zelig(unem ~ gdp + capmob + trade, model = "normal.bayes",
+ data = macro, verbose=TRUE)
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+
+heidel.diag(z.out$coefficients)
+user.prompt()
+
+raftery.diag(z.out$coefficients)
+user.prompt()
+
+## summarizing the output
+summary(z.out)
+user.prompt()
+
+## Setting values for the explanatory variables to
+## their sample averages:
+x.out <- setx(z.out)
+user.prompt()
+
+## Simulating quantities of interest from the posterior
+## distribution given x.out:
+s.out1 <- sim(z.out, x = x.out)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out1)
+user.prompt()
+
+## Simulating First Differences
+## Set explanatory variables to their default(mean/mode) values,
+## with high (80th percentile) and low (20th percentile) trade on GDP:
+x.high <- setx(z.out, trade = quantile(macro$trade, prob = 0.8))
+x.low <- setx(z.out, trade = quantile(macro$trade, prob = 0.2))
+user.prompt()
+
+## Estimating the first difference for the effect of
+## high versus low trade on unemployment rate:
+s.out2 <- sim(z.out, x = x.high, x1 = x.low)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out2)
+
+
+
+
+
+
+
+
diff --git a/demo/oprobit.bayes.R b/demo/oprobit.bayes.R
new file mode 100644
index 0000000..7682503
--- /dev/null
+++ b/demo/oprobit.bayes.R
@@ -0,0 +1,67 @@
+## Attaching the example dataset:
+data(sanction)
+
+# Create an ordered dependent variable:
+user.prompt()
+sanction$ncost <- factor(sanction$ncost, ordered = TRUE,
+ levels = c("net gain", "little effect",
+ "modest loss", "major loss"))
+
+## Estimating the model using oprobit.bayes:
+z.out <- zelig(ncost ~ mil + coop, model = "oprobit.bayes",
+ data = sanction, verbose=TRUE)
+
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+#geweke.diag(z.out$coefficients)
+#user.prompt()
+
+heidel.diag(z.out$coefficients)
+user.prompt()
+
+raftery.diag(z.out$coefficients)
+user.prompt()
+
+## summarizing the output
+summary(z.out)
+user.prompt()
+
+## Setting values for the explanatory variables to
+## their sample averages:
+x.out <- setx(z.out)
+user.prompt()
+
+## Simulating quantities of interest from the posterior
+## distribution given x.out:
+s.out1 <- sim(z.out, x = x.out)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out1)
+user.prompt()
+
+## Simulating First Differences:
+## Setting explanatory variables to their default(mean/mode)
+## values, with military action to be yes(1) or no(0)
+x.high <- setx(z.out, mil=0)
+x.low <- setx(z.out, mil=1)
+user.prompt()
+
+## Estimating the first difference for the effect of
+## military action on the probabilities of
+## incurring differnt level of cost:
+
+s.out2 <- sim(z.out, x = x.high, x1 = x.low)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out2)
+
+
+
+
+
+
+
+
diff --git a/demo/poisson.bayes.R b/demo/poisson.bayes.R
new file mode 100644
index 0000000..81b3db6
--- /dev/null
+++ b/demo/poisson.bayes.R
@@ -0,0 +1,59 @@
+## Attaching the example dataset:
+data(sanction)
+
+## Estimating the model using poisson.bayes:
+z.out <- zelig(num ~ target + coop, model = "poisson.bayes",
+ data = sanction, verbose=TRUE)
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+
+heidel.diag(z.out$coefficients)
+user.prompt()
+
+raftery.diag(z.out$coefficients)
+user.prompt()
+
+## summarizing the output
+summary(z.out)
+user.prompt()
+
+## Setting values for the explanatory variables to
+## their sample averages:
+x.out <- setx(z.out)
+user.prompt()
+
+## Simulating quantities of interest from the posterior
+## distribution given x.out:
+s.out1 <- sim(z.out, x = x.out)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out1)
+user.prompt()
+
+## Simulating First Differences:
+## Setting explanatory variables to their default(mean/mode)
+## values, with the number of targets to be its maximum
+## versus its minimum:
+x.max <- setx(z.out, target = max(sanction$target))
+x.min <- setx(z.out, target = min(sanction$target))
+user.prompt()
+
+## Estimating the first difference for the effect of
+## maximum versus minimum number of targets:
+s.out2 <- sim(z.out, x = x.max, x1 = x.min)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out2)
+
+
+
+
+
+
+
+
diff --git a/demo/probit.bayes.R b/demo/probit.bayes.R
new file mode 100644
index 0000000..40b87ad
--- /dev/null
+++ b/demo/probit.bayes.R
@@ -0,0 +1,59 @@
+## Attaching the example dataset:
+data(turnout)
+
+## Estimating the model using probit.bayes:
+z.out <- zelig(vote ~ race + educate, model = "probit.bayes",
+ data = turnout, verbose=TRUE)
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+
+heidel.diag(z.out$coefficients)
+user.prompt()
+
+raftery.diag(z.out$coefficients)
+user.prompt()
+
+## summarizing the output
+summary(z.out)
+user.prompt()
+
+## Setting values for the explanatory variables to
+## their sample averages:
+x.out <- setx(z.out)
+user.prompt()
+
+## Simulating quantities of interest from the posterior
+## distribution given x.out:
+s.out1 <- sim(z.out, x = x.out)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out1)
+user.prompt()
+
+## Simulating First Differences:
+## Setting education is set to be between low(25th percentile)
+## versus high(75th percentile) while all the other variables
+## held at their default values.
+x.high <- setx(z.out, educate = quantile(turnout$educate, prob = 0.75))
+x.low <- setx(z.out, educate = quantile(turnout$educate, prob = 0.25))
+user.prompt()
+
+## Estimating the first difference for the effect of
+## high versus low trade on unemployment rate:
+s.out2 <- sim(z.out, x = x.high, x1 = x.low)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out2)
+
+
+
+
+
+
+
+
diff --git a/demo/relogit.R b/demo/relogit.R
index 6f41180..894b226 100644
--- a/demo/relogit.R
+++ b/demo/relogit.R
@@ -1,28 +1,59 @@
data(mid)
user.prompt()
+
+
+## prior correction + bias correction
z.out1 <- zelig(conflict ~ major + contig + power + maxdem + mindem + years,
data = mid, model = "relogit", tau = 1042/303772)
user.prompt()
+
summary(z.out1)
user.prompt()
+
x.out1 <- setx(z.out1)
user.prompt()
+
s.out1 <- sim(z.out1, x = x.out1)
user.prompt()
+
summary(s.out1)
user.prompt()
+
plot(s.out1)
-user.prompt()
+
+## weighting + bias correction + robust s.e.
z.out2 <- zelig(conflict ~ major + contig + power + maxdem + mindem + years,
- data = mid, model = "relogit", tau = c(0.002, 0.005))
+ data = mid, model = "relogit", tau = 1042/303772,
+ case.control = "weighting", robust = TRUE)
user.prompt()
+
summary(z.out2)
user.prompt()
+
x.out2 <- setx(z.out2)
user.prompt()
+
s.out2 <- sim(z.out2, x = x.out2)
user.prompt()
-summary(s.out2)
+
+
+## bounds
+user.prompt()
+z.out3 <- zelig(conflict ~ major + contig + power + maxdem + mindem + years,
+ data = mid, model = "relogit", tau = c(0.002, 0.005))
user.prompt()
-plot(s.out2)
+
+summary(z.out3)
+user.prompt()
+
+x.out3 <- setx(z.out3)
+user.prompt()
+
+s.out3 <- sim(z.out3, x = x.out3)
+user.prompt()
+
+summary(s.out3)
+user.prompt()
+
+plot(s.out3)
diff --git a/demo/tobit.bayes.R b/demo/tobit.bayes.R
new file mode 100644
index 0000000..30726a3
--- /dev/null
+++ b/demo/tobit.bayes.R
@@ -0,0 +1,59 @@
+## Attaching the example dataset:
+data(tobin)
+
+## Estimating the model using tobit.bayes:
+z.out <- zelig(durable ~ age + quant, model = "tobit.bayes",
+ data = tobin, verbose=TRUE)
+user.prompt()
+
+## Checking for convergence before summarizing the estimates:
+geweke.diag(z.out$coefficients)
+user.prompt()
+
+heidel.diag(z.out$coefficients)
+user.prompt()
+
+raftery.diag(z.out$coefficients)
+user.prompt()
+
+## summarizing the output
+summary(z.out)
+user.prompt()
+
+## Setting values for the explanatory variables to
+## their sample averages:
+x.out <- setx(z.out)
+user.prompt()
+
+## Simulating quantities of interest from the posterior
+## distribution given x.out:
+s.out1 <- sim(z.out, x = x.out)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out1)
+user.prompt()
+
+## Simulating First Differences:
+## Setting explanatory variables to their default(mean/mode)
+## values, with high (80th percentile) and low (20th percentile)
+## liquidity ratio(\texttt{quant}):
+x.high <- setx(z.out, quant = quantile(tobin$quant, prob = 0.8))
+x.low <- setx(z.out, quant = quantile(tobin$quant, prob = 0.2))
+user.prompt()
+
+## Estimating the first difference for the effect of
+## high versus low liquidity ratio:
+s.out2 <- sim(z.out, x = x.high, x1 = x.low)
+user.prompt()
+
+## Summarizing the simulation results:
+summary(s.out2)
+
+
+
+
+
+
+
+
diff --git a/man/PErisk.Rd b/man/PErisk.Rd
new file mode 100644
index 0000000..1655dc3
--- /dev/null
+++ b/man/PErisk.Rd
@@ -0,0 +1,76 @@
+\name{PErisk}
+
+\alias{PErisk}
+
+\title{Political Economic Risk Data from 62 Countries in 1987}
+
+\description{
+ Political Economic Risk Data from 62 Countries in 1987.
+
+}
+
+\usage{data(PErisk)}
+
+\format{
+ A data frame with 62 observations on the following 6 variables.
+ All data points are from 1987. See Quinn (2004) for more
+ details.
+
+ country: a factor with levels 'Argentina' 'Australia' 'Austria'
+ 'Bangladesh' 'Belgium' 'Bolivia' 'Botswana' 'Brazil' 'Burma'
+ 'Cameroon' 'Canada' 'Chile' 'Colombia' 'Congo-Kinshasa'
+ 'Costa Rica' 'Cote d'Ivoire' 'Denmark' 'Dominican Republic'
+ 'Ecuador' 'Finland' 'Gambia, The' 'Ghana' 'Greece' 'Hungary'
+ 'India' 'Indonesia' 'Iran' 'Ireland' 'Israel' 'Italy' 'Japan'
+ 'Kenya' 'Korea, South' 'Malawi' 'Malaysia' 'Mexico' 'Morocco'
+ 'New Zealand' 'Nigeria' 'Norway' 'Papua New Guinea'
+ 'Paraguay' 'Philippines' 'Poland' 'Portugal' 'Sierra Leone'
+ 'Singapore' 'South Africa' 'Spain' 'Sri Lanka' 'Sweden'
+ 'Switzerland' 'Syria' 'Thailand' 'Togo' 'Tunisia' 'Turkey'
+ 'United Kingdom' 'Uruguay' 'Venezuela' 'Zambia' 'Zimbabwe'
+
+ courts: an ordered factor with levels '0' < '1'.'courts' is an
+ indicator of whether the country in question is judged to
+ have an independent judiciary. From Henisz (2002).
+
+ barb2: a numeric vector giving the natural log of the black market
+ premium in each country. The black market premium is coded as
+ the black market exchange rate (local currency per dollar)
+ divided by the official exchange rate minus 1. From
+ Marshall, Gurr, and Harff (2002).
+
+ prsexp2: an ordered factor with levels '0' < '1' < '2' < '3' < '4'
+ < '5', giving the lack of expropriation risk. From Marshall,
+ Gurr, and Harff (2002).
+
+ prscorr2: an ordered factor with levels '0' < '1' < '2' < '3' < '4'
+ < '5', measuring the lack of corruption. From Marshall, Gurr,
+ and Harff (2002).
+
+ gdpw2: a numeric vector giving the natural log of real GDP per
+ worker in 1985 international prices. From Alvarez et al.
+ (1999).
+}
+
+\source{
+ Mike Alvarez, Jose Antonio Cheibub, Fernando Limongi, and Adam
+ Przeworski. 1999. ``ACLP Political and Economic Database.'' <URL:
+ http://www.ssc.upenn.edu/~cheibub/data/>.
+
+ Witold J. Henisz. 2002. ``The Political Constraint Index (POLCON)
+ Dataset.'' \ <URL:
+ http://www-management.wharton.upenn.edu/henisz/POLCON/ContactInfo.
+ html>.
+
+ Monty G. Marshall, Ted Robert Gurr, and Barbara Harff. 2002.
+ ``State Failure Task Force Problem Set.'' <URL:
+ http://www.cidcm.umd.edu/inscr/stfail/index.htm>.
+}
+
+\references{
+ Kevin M. Quinn. 2004. ``Bayesian Factor Analysis for Mixed Ordinal
+ and Continuous Response.'' \emph{Political Analyis}. Vol. 12, pp.338--353.
+}
+
+
+\keyword{datasets}
diff --git a/man/SupremeCourt.Rd b/man/SupremeCourt.Rd
new file mode 100644
index 0000000..9b04226
--- /dev/null
+++ b/man/SupremeCourt.Rd
@@ -0,0 +1,31 @@
+\name{SupremeCourt}
+
+\alias{SupremeCourt}
+
+\title{U.S. Supreme Court Vote Matrix}
+
+\description{
+ This dataframe contains a matrix votes cast by U.S. Supreme Court
+ justices in all cases in the 2000 term.
+}
+
+\usage{data(SupremeCourt)}
+
+\format{ The dataframe has contains data for justices Rehnquist, Stevens,
+ O'Connor, Scalia, Kennedy, Souter, Thomas, Ginsburg, and Breyer
+ for the 2000 term of the U.S. Supreme Court. It contains data
+ from 43 non-unanimous cases. The votes are coded liberal (1) and
+ conservative (0) using the protocol of Spaeth (2003). The unit
+ of analysis is the case citation (ANALU=0). We are concerned with
+ formally decided cases issued with written opinions, after full
+ oral argument and cases decided by an equally divided vote
+ (DECTYPE=1,5,6,7).}
+
+\source{
+ Harold J. Spaeth (2005). ``Original United States Supreme Court
+ Database: 1953-2004 Terms.''
+ <URL:http://www.as.uky.edu/polisci/ulmerproject/sctdata.htm>.
+}
+
+
+\keyword{datasets}
diff --git a/man/eidat.Rd b/man/eidat.Rd
new file mode 100644
index 0000000..c953306
--- /dev/null
+++ b/man/eidat.Rd
@@ -0,0 +1,19 @@
+\name{eidat}
+
+\alias{eidat}
+
+\title{Simulation Data for Ecological Inference}
+
+\description{
+ This dataframe contains a simulated data set to illustrate the models
+ for ecological inference.
+}
+
+\usage{data(eidat)}
+
+\format{
+ A table containing 4 variables ("t0", "t1", "x0", "x1") and 10
+observations.
+}
+
+\keyword{datasets}
diff --git a/man/newpainters.Rd b/man/newpainters.Rd
new file mode 100644
index 0000000..0995248
--- /dev/null
+++ b/man/newpainters.Rd
@@ -0,0 +1,40 @@
+\name{newpainters}
+
+\alias{newpainters}
+
+\title{The Discretized Painter's Data of de Piles}
+
+\description{
+ The original painters data contain the subjective assessment,
+ on a 0 to 20 integer scale, of 54 classical painters. The
+ newpainters data discretizes the subjective assessment by
+ quartiles with thresholds 25\%, 50\%, 75\%. The painters were
+ assessed on four characteristics: composition, drawing,
+ colour and expression. The data is due to the Eighteenth century
+ art critic, de Piles.
+
+}
+
+\usage{data(newpainters)}
+
+\format{A table containing 5 variables ("Composition", "Drawing", "Colour",
+"Expression", and "School") and 54 observations.}
+
+\source{
+
+ A. J. Weekes (1986).``A Genstat Primer''. Edward Arnold.
+
+ M. Davenport and G. Studdert-Kennedy (1972). ``The statistical
+ analysis of aesthetic judgement: an exploration.'' \emph{Applied
+ Statistics}, vol. 21, pp. 324--333.
+
+ I. T. Jolliffe (1986) ``Principal Component Analysis.'' Springer.
+}
+
+\references{
+
+ Venables, W. N. and Ripley, B. D. (2002) ``Modern Applied
+ Statistics with S,'' Fourth edition. Springer.
+}
+
+\keyword{datasets}
diff --git a/man/swiss.Rd b/man/swiss.Rd
new file mode 100644
index 0000000..ab88b9a
--- /dev/null
+++ b/man/swiss.Rd
@@ -0,0 +1,51 @@
+\name{swiss}
+
+\alias{swiss}
+
+\title{Swiss Fertility and Socioeconomic Indicators (1888) Data}
+
+\description{
+ Standardized fertility measure and socio-economic indicators for
+ each of 47 French-speaking provinces of Switzerland at about 1888.
+}
+
+\usage{data(swiss)}
+
+\format{
+ A data frame with 47 observations on 6 variables, each of which
+ is in percent, i.e., in [0,100].
+
+ [,1] Fertility Ig, "common standardized fertility measure"
+ [,2] Agriculture % of males involved in agriculture as occupation
+ [,3] Examination % "draftees" receiving highest mark on army exami
+nation
+ [,4] Education % education beyond primary school for "draftees".
+ [,5] Catholic % catholic (as opposed to "protestant").
+ [,6] Infant.Mortality live births who live less than 1 year.
+
+ All variables but 'Fert' give proportions of the population.
+}
+
+\source{
+
+
+ Project "16P5", pages 549-551 in
+
+ Mosteller, F. and Tukey, J. W. (1977) ``Data Analysis and
+ Regression: A Second Course in Statistics''. Addison-Wesley,
+ Reading Mass.
+
+ indicating their source as "Data used by permission of Franice van
+ de Walle. Office of Population Research, Princeton University,
+ 1976. Unpublished data assembled under NICHD contract number No
+ 1-HD-O-2077."
+ }
+
+ \references{
+ Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) ``The New S
+ Language''. Wadsworth & Brooks/Cole.
+
+}
+
+
+\keyword{datasets}
diff --git a/man/tobin.Rd b/man/tobin.Rd
new file mode 100644
index 0000000..3f40806
--- /dev/null
+++ b/man/tobin.Rd
@@ -0,0 +1,32 @@
+\name{tobin}
+
+\alias{tobin}
+
+\title{Tobin's Tobit Data}
+
+\description{
+ Economists fit a parametric censored data model called the
+ `tobit'. These data are from Tobin's original paper.
+}
+
+\usage{data(tobin)}
+
+\format{
+ A data frame with 20 observations on the following 3 variables.
+
+ durable: Durable goods purchase
+
+ age: Age in years
+
+ quant: Liquidity ratio (x 1000)
+}
+
+\source{
+ J. Tobin, Estimation of relationships for limited dependent
+ variables, Econometrica, v26, 24-36, 1958.
+ }
+
+}
+
+
+\keyword{datasets}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-zelig.git
More information about the debian-science-commits
mailing list