[r-cran-zelig] 08/102: Import Upstream version 2.2-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 1a5deb70b76bff79cc828369ee67425c3eb7e975
Author: Andreas Tille <tille at debian.org>
Date: Sun Jan 8 09:38:58 2017 +0100
Import Upstream version 2.2-1
---
DESCRIPTION | 9 +++---
NAMESPACE | 33 +++++++++++++-------
R/bootfn.default.R | 2 +-
R/define.data.R | 34 ---------------------
R/define.par.R | 20 ------------
R/formula.vglm.R | 3 ++
R/generate.start.R | 9 ++++++
R/model.end.R | 29 +++++++-----------
R/param.vglm.R | 2 +-
R/parse.par.R | 20 +++---------
R/print.setx.R | 7 -----
R/print.summary.glm.robust.R | 71 +++++++++++++++++++++++++++++++++++++++++++
R/print.summary.lm.robust.R | 72 ++++++++++++++++++++++++++++++++++++++++++++
R/print.summary.zelig.R | 8 ++---
R/qi.polr.R | 52 +++++++++++++++-----------------
R/qi.vglm.R | 46 ++++++----------------------
R/relogit.R | 4 +--
R/setx.MI.R | 4 +--
R/setx.default.R | 48 ++++++++++++-----------------
R/setx.strata.R | 2 +-
R/sim.cond.R | 28 +++++------------
R/sim.default.R | 3 ++
R/summary.glm.robust.R | 6 +++-
R/summary.lm.robust.R | 11 +++++--
R/zelig.R | 40 ++++++------------------
R/zelig2beta.R | 2 +-
R/zelig2blogit.R | 3 +-
R/zelig2bprobit.R | 2 +-
R/zelig2exp.R | 21 +++++++++++--
R/zelig2gamma.R | 6 ++--
R/zelig2logit.R | 6 ++--
R/zelig2lognorm.R | 21 +++++++++++--
R/zelig2ls.R | 4 +--
R/zelig2mlogit.R | 4 +--
R/zelig2mloglm.R | 2 +-
R/zelig2negbin.R | 2 +-
R/zelig2normal.R | 6 ++--
R/zelig2ologit.R | 3 +-
R/zelig2oprobit.R | 27 +++--------------
R/zelig2poisson.R | 6 ++--
R/zelig2probit.R | 7 ++---
R/zelig2weibull.R | 21 +++++++++++--
R/zelig3glm.R | 22 ++++++++++++++
R/zelig3ls.R | 19 ++++++++++++
R/zelig3ologit.R | 12 ++++++++
R/zelig3oprobit.R | 12 ++++++++
README | 10 ++++--
demo/00Index | 1 +
demo/oprobit.R | 9 +++---
demo/robust.R | 33 ++++++++++++++++++--
demo/vertci.R | 6 ++--
man/help.zelig.Rd | 14 +++------
man/summary.zelig.Rd | 2 +-
53 files changed, 500 insertions(+), 346 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 050394a..9919bbf 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,12 +1,13 @@
Package: Zelig
-Version: 2.1-4
-Date: 2005-05-22
+Version: 2.2-1
+Date: 2005-7-11
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
+Depends: R (>= 1.9.1), MASS, boot,
+Suggests: VGAM, survival, sandwich
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
@@ -22,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: Sun May 22 13:46:26 2005; king
+Packaged: Mon Jul 11 12:39:01 2005; olau
diff --git a/NAMESPACE b/NAMESPACE
index e3707fe..d70f148 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,24 +1,34 @@
-import(MASS)
-import(boot)
+importFrom(MASS, mvrnorm)
+importFrom(boot, boot)
+importFrom(stats, glm)
+importFrom(stats, family)
+importFrom(stats, binomial)
+importFrom(stats, gaussian)
+importFrom(stats, Gamma)
+importFrom(stats, inverse.gaussian)
+importFrom(stats, poisson)
+importFrom(stats, quasi)
+importFrom(stats, quasibinomial)
+importFrom(stats, quasipoisson)
export( dims,
- help.zelig,
- zelig,
- help.zelig,
- gsource,
+ gsource,
+ help.zelig,
+ plot.ci,
+ plot.zelig,
+ rocplot,
setx,
sim,
- plot.ci,
- plot.zelig,
- rocplot,
ternarypoints,
- user.prompt
+ user.prompt,
+ zelig
)
S3method("$", vglm)
S3method("$", summary.vglm)
S3method("$<-", vglm)
S3method("$<-", summary.vglm)
+S3method(formula, vglm)
S3method(names, relogit)
S3method(names, summary.zelig.relogit)
S3method(names, vglm)
@@ -46,10 +56,11 @@ S3method(print, BetaReg)
S3method(print, names.relogit)
S3method(print, names.zelig)
S3method(print, relogit)
-S3method(print, setx)
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.zelig)
S3method(print, zelig)
S3method(repl, zelig)
diff --git a/R/bootfn.default.R b/R/bootfn.default.R
index 32faa27..6dcbef3 100644
--- a/R/bootfn.default.R
+++ b/R/bootfn.default.R
@@ -1,8 +1,8 @@
bootfn.default <- function(data, i, object) {
d <- data[i,]
object$call$data <- d
- fit <- eval(object$call, sys.parent())
l <- length(param(object, bootstrap = TRUE))
+ fit <- eval(object$call, sys.parent())
l1 <- length(param(fit, bootstrap = TRUE))
while (l > l1) {
object$call$data <- data[sample(nrow(data), replace=TRUE),]
diff --git a/R/define.data.R b/R/define.data.R
deleted file mode 100644
index 50245ce..0000000
--- a/R/define.data.R
+++ /dev/null
@@ -1,34 +0,0 @@
-define.data <- function(formula, data, x.name = "x",
- y.name = "y", mf.name = "mf",
- xlev.name = "xlev", lev.name = "lev") {
- mc <- match.call()
- mc$x.name <- mc$y.name <- mc$t.name <- mc$xlev.name <- mc$l.name <- NULL
-
- mc[[1]] <- as.name("model.frame")
- mf <- eval.parent(mc)
- assign(mf.name, mf, envir = parent.frame())
-
- Terms <- attr(mf, "terms")
-
- x <- model.matrix(Terms, mf)
- assign(x.name, x, envir = parent.frame())
-
- yvar <- attr(Terms, "response")
- xvars <- as.character(attr(Terms, "response"))[-1]
- xlev <- if (length(xvars) > 0) {
- xlev <- lapply(mf[xvars], levels)
- xlev[!sapply(xlev, is.null)]
- }
- assign(xlev.name, xlev, envir = parent.frame())
-
- y <- model.response(mf)
- assign(y.name, y, envir = parent.frame())
-
- lev <- levels(y)
- assign(lev.name, lev, envir = parent.frame())
-
- attr(Terms, "items") <- list(x = x.name, y = y.name,
- mf = mf.name, xlev =
- xlev.name, lev = lev.name)
- Terms
-}
diff --git a/R/define.par.R b/R/define.par.R
deleted file mode 100644
index 101a43d..0000000
--- a/R/define.par.R
+++ /dev/null
@@ -1,20 +0,0 @@
-define.par <- function(object, coef.name = "beta",
- ancillary = NULL, start = NULL){
- mc <- match.call()
- x <- eval(as.name(attr(object, "items")$x), parent.frame())
- coef <- colnames(x)
- if (!is.null(ancillary))
- coef <- c(coef, ancillary)
-
- attr(object, "items")$ancillary <- ancillary
- attr(object, "items")$coef <- coef.name
- assign(as.character(mc$object), object, envir = parent.frame())
-
- if (is.null(start))
- start.val <- vector(mode = "numeric", length = length(coef))
- else
- start.val <- start
-
- names(start.val) <- coef
- start.val
-}
diff --git a/R/formula.vglm.R b/R/formula.vglm.R
new file mode 100644
index 0000000..4d61970
--- /dev/null
+++ b/R/formula.vglm.R
@@ -0,0 +1,3 @@
+formula.vglm <- function(x, ...) {
+ x at call$formula
+}
diff --git a/R/generate.start.R b/R/generate.start.R
new file mode 100644
index 0000000..b509bf0
--- /dev/null
+++ b/R/generate.start.R
@@ -0,0 +1,9 @@
+generate.start <- function(start.val, X, ancillary = NULL) {
+ if (is.null(start.val))
+ start.val <- rep(0, ncol(X) + length(ancillary))
+ if (!is.null(ancillary))
+ names(start.val) <- c(colnames(X), ancillary)
+ else
+ names(start.val) <- c(colnames(X))
+ start.val
+}
diff --git a/R/model.end.R b/R/model.end.R
index fb38309..c02e42f 100644
--- a/R/model.end.R
+++ b/R/model.end.R
@@ -1,24 +1,17 @@
-model.end <- function(Terms, res) {
+model.end <- function(res, D) {
- variance <- -solve(res$hessian)
+ res$variance <- -solve(res$hessian)
+ res$hessian <- NULL
- lev <- eval(as.name(attr(Terms, "items")$lev), envir = parent.frame())
- kall <- eval(as.name("kall"), envir = parent.frame())
- xlev <- eval(as.name(attr(Terms, "items")$xlev), envir = parent.frame())
- mf <- eval(as.name(attr(Terms, "items")$mf), envir = parent.frame())
- x <- eval(as.name(attr(Terms, "items")$x), envir = parent.frame())
-
- names(res$par) <- c(colnames(x), attr(Terms, "items")$ancillary)
- colnames(variance) <- rownames(variance) <- names(res$par)
+ colnames(res$variance) <- rownames(res$variance) <- names(res$par)
+ res$coefficients <- res$par
+ res$par <- NULL
- fit <- list(coefficients = res$par, variance = variance,
- lev = lev, terms = Terms, call = kall,
- convergence = res$convergence, xlevels = xlev)
+ res$terms <- attr(D, "terms")
- attr(fit, "na.message") <- attr(mf, "na.message")
+ attr(res, "na.message") <- attr(D, "na.message")
+ if (!is.null(attr(D, "na.action")))
+ res$na.action <- attr(D, "na.action")
- if (!is.null(attr(mf, "na.action")))
- fit$na.action <- attr(mf, "na.action")
-
- fit
+ res
}
diff --git a/R/param.vglm.R b/R/param.vglm.R
index 99c987b..ccbace9 100644
--- a/R/param.vglm.R
+++ b/R/param.vglm.R
@@ -1,5 +1,5 @@
param.vglm <- function(object, num, bootstrap = FALSE) {
- cov <- vcov(object)
+ cov <- VGAM::vcov(object)
res <- object at coefficients
if (!bootstrap)
res <- mvrnorm(num, mu=res, Sigma=cov)
diff --git a/R/parse.par.R b/R/parse.par.R
index 68ae68d..15dea51 100644
--- a/R/parse.par.R
+++ b/R/parse.par.R
@@ -1,17 +1,7 @@
-parse.par <- function(par, x, terms) {
- k <- ncol(x)
- beta <- par[1:k]
- assign(attr(terms, "items")$coef, beta, envir = parent.frame())
- ancil.names <- attr(terms, "items")$ancillary
- j <- 0
- if (!is.null(ancil.names)) {
- j <- length(ancil.names)
- for (i in 1:j) {
- tmp <- par[(k + i)]
- assign(ancil.names[i], tmp, envir = parent.frame())
- }
- }
- if (length(par) > (k + j))
- return(par[(k + j):length(par)])
+parse.par <- function(par, idx, name) {
+ if (name %in% idx)
+ return(par[which(idx == name)])
+ else
+ return(par[which(idx != name)])
}
diff --git a/R/print.setx.R b/R/print.setx.R
deleted file mode 100644
index 9fd17bb..0000000
--- a/R/print.setx.R
+++ /dev/null
@@ -1,7 +0,0 @@
-print.setx<-function(x, digits=getOption("digits"),
- ...){
- attributes(x)$class<-NULL
- attributes(x)$assign<-NULL
- attributes(x)$contrasts<-NULL
- print.matrix(x, digits=digits, ...)
-}
diff --git a/R/print.summary.glm.robust.R b/R/print.summary.glm.robust.R
new file mode 100644
index 0000000..e26134c
--- /dev/null
+++ b/R/print.summary.glm.robust.R
@@ -0,0 +1,71 @@
+print.summary.glm.robust <-
+ function (x, digits = max(3, getOption("digits") - 3),
+ 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)
+}
+
diff --git a/R/print.summary.lm.robust.R b/R/print.summary.lm.robust.R
new file mode 100644
index 0000000..f8aa465
--- /dev/null
+++ b/R/print.summary.lm.robust.R
@@ -0,0 +1,72 @@
+print.summary.lm.robust <-
+ function (x, digits = max(3, getOption("digits") - 3),
+ symbolic.cor = x$symbolic.cor,
+ signif.stars= getOption("show.signif.stars"), ...)
+{
+ cat("\nCall:\n")#S: ' ' instead of '\n'
+ cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")
+ resid <- x$residuals
+ df <- x$df
+ rdf <- df[2]
+ cat(if(!is.null(x$w) && diff(range(x$w))) "Weighted ",
+ "Residuals:\n", sep="")
+ if (rdf > 5) {
+ nam <- c("Min", "1Q", "Median", "3Q", "Max")
+ rq <- if (length(dim(resid)) == 2)
+ structure(apply(t(resid), 1, quantile),
+ dimnames = list(nam, dimnames(resid)[[2]]))
+ else structure(quantile(resid), names = nam)
+ print(rq, digits = digits, ...)
+ }
+ else if (rdf > 0) {
+ print(resid, digits = digits, ...)
+ } else { # rdf == 0 : perfect fit!
+ cat("ALL", df[1], "residuals are 0: no residual degrees of freedom!\n")
+ }
+ if (length(x$aliased) == 0) {
+ cat("\nNo Coefficients\n")
+ } else {
+ if (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("\nResidual standard error:",
+ format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")
+ if (!is.null(x$fstatistic)) {
+ cat("Multiple R-Squared:", formatC(x$r.squared, digits=digits))
+ cat(",\tAdjusted R-squared:",formatC(x$adj.r.squared,digits=digits),
+ "\nF-statistic:", formatC(x$fstatistic[1], digits=digits),
+ "on", x$fstatistic[2], "and",
+ x$fstatistic[3], "DF, p-value:",
+ format.pval(pf(x$fstatistic[1], x$fstatistic[2],
+ x$fstatistic[3], lower.tail = FALSE), digits=digits),
+ "\n")
+ }
+ correl <- x$correlation
+ if (!is.null(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")#- not in S
+ invisible(x)
+}
diff --git a/R/print.summary.zelig.R b/R/print.summary.zelig.R
index b081acb..629cde8 100644
--- a/R/print.summary.zelig.R
+++ b/R/print.summary.zelig.R
@@ -9,10 +9,10 @@ print.summary.zelig <- function(x, digits=getOption("digits"),
cat("\nObserved Data \n")
else
cat("\nValues of X \n")
- print.setx(x$x, digits=digits, ...)
+ print(x$x, digits=digits, ...)
if(!is.null(x$x1)){
cat("\nValues of X1 \n")
- print.setx(x$x1, digits=digits, ...)
+ print(x$x1, digits=digits, ...)
}
}
else {
@@ -20,10 +20,10 @@ print.summary.zelig <- function(x, digits=getOption("digits"),
cat("\nMean Values of Observed Data (n = ", nrow(x$x), ") \n", sep = "")
else
cat("\nMean Values of X (n = ", nrow(x$x), ") \n", sep = "")
- print.setx(apply(x$x, 2, mean), digits=digits, ...)
+ print(apply(x$x, 2, mean), digits=digits, ...)
if (!is.null(x$x1)) {
cat("\nMean Values of X1 (n = ", nrow(x$x1), ") \n", sep = "")
- print.setx(apply(x$x1, 2, mean), digits=digits, ...)
+ print(apply(x$x1, 2, mean), digits=digits, ...)
}
}
}
diff --git a/R/qi.polr.R b/R/qi.polr.R
index 9467f71..d6a13e9 100644
--- a/R/qi.polr.R
+++ b/R/qi.polr.R
@@ -1,42 +1,38 @@
qi.polr <- function(object, simpar, x, x1 = NULL, y = NULL) {
+ num <- nrow(simpar)
m <- length(coef(object))
sim.coef <- simpar[,1:m]
sim.zeta <- simpar[,(m+1):ncol(simpar)]
k <- length(object$zeta) + 1
lev <- object$lev
- eta <- x[,-1] %*% t(sim.coef)
- ev.polr <- function(ev, eta, sim.zeta) {
- for (j in 1:ncol(sim.zeta))
- ev[,j,] <- exp(sim.zeta[,j] - t(eta)) / (1 + exp(sim.zeta[,j] - t(eta)))
- for (j in 0:(ncol(sim.zeta)-1))
- ev[,dim(ev)[2]-j,] <- ev[,dim(ev)[2]-j,]-ev[,(dim(ev)[2]-j-1),]
- ev
- }
- ev <- array(1, dim = c(nrow(sim.coef), length(lev), nrow(x)),
- dimnames = list(NULL, lev, rownames(x)))
- ev <- ev.polr(ev, eta, sim.zeta)
- sim.cut <- Ipr <- array(1, dim = dim(ev), dimnames = dimnames(ev))
- pr <- matrix(NA, nrow = nrow(sim.coef), ncol = nrow(x))
- sim.cut[,1,] <- ev[,1,]
- tmp <- matrix(runif(length(sim.cut[,1,]), 0, 1), nrow =
- dim(sim.cut)[1], ncol = dim(sim.cut)[3])
- for (j in 2:length(lev))
- sim.cut[,j,] <- sim.cut[,(j-1),] + ev[,j,]
- for (k in 1:length(lev))
- Ipr[,k,] <- as.integer(tmp > sim.cut[,k,])
- for (n in 1:dim(Ipr)[3])
+ eta <- t(x[,-1] %*% t(sim.coef))
+ Ipr <- cuts <- tmp0 <- array(0, dim = c(num, k, nrow(x)),
+ dimnames = list(1:num, lev, rownames(x)))
+ for (i in 1:num)
+ cuts[i,,] <- t(object$inv.link(eta[i,], sim.zeta[i,]))
+ tmp0[,(2:k),] <- cuts[,(1:(k-1)),]
+ ev <- cuts - tmp0
+ pr <- matrix(NA, nrow = num, ncol = nrow(x))
+ tmp <- matrix(runif(length(cuts[,1,]), 0, 1),
+ nrow = num , ncol = nrow(x))
+ for (i in 1:k)
+ Ipr[,i,] <- as.integer(tmp > cuts[,i,])
+ for (n in 1:nrow(x))
pr[,n] <- 1 + rowSums(Ipr[,,n])
pr <- matrix(factor(pr, labels = lev, ordered = TRUE),
nrow = nrow(sim.coef), ncol = nrow(x))
colnames(pr) <- rownames(x)
- qi <- list(ev=ev, pr=pr)
- qi.name <- list(ev="Expected Values: P(Y=j|X)",
- pr="Predicted Values: Y|X")
+ qi <- list(ev = ev, pr = pr)
+ qi.name <- list(ev = "Expected Values: P(Y=j|X)",
+ pr = "Predicted Values: Y|X")
if(!is.null(x1)){
- ev1 <- array(1, dim = c(nrow(sim.coef), length(lev), nrow(x)),
- dimnames = list(NULL, lev, rownames(x)))
- eta1 <- x1[,-1] %*% t(sim.coef)
- ev1 <- ev.polr(ev1, eta1, sim.zeta)
+ eta1 <- t(x1[,-1] %*% t(sim.coef))
+ Ipr <- cuts <- tmp0 <- array(0, dim = c(num, k, nrow(x)),
+ dimnames = list(1:num, lev, rownames(x)))
+ for (i in 1:num)
+ cuts[i,,] <- t(object$inv.link(eta1[i,], sim.zeta[i,]))
+ tmp0[,(2:k),] <- cuts[,(1:(k-1)),]
+ ev1 <- cuts - tmp0
qi$fd <- ev1 - ev
qi$rr <- ev1 / ev
qi.name$fd <- "First Differences: P(Y=j|X1)-P(Y=j|X)"
diff --git a/R/qi.vglm.R b/R/qi.vglm.R
index f9297ba..fe223f8 100644
--- a/R/qi.vglm.R
+++ b/R/qi.vglm.R
@@ -1,12 +1,10 @@
qi.vglm <- function (object, simpar, x, x1=NULL, y = NULL) {
model <- object$zelig
cm <- object at constraints
- if (model=="mlogit" || model=="oprobit")
+ if (model=="mlogit")
ndim <- (ncol(object at y)-1)
else if (model=="blogit" || model=="bprobit")
ndim <- 3
- else if (model=="sur")
- ndim <- ncol(object at y)
v <- rep(list(NULL), ndim)
for(i in 1:length(cm)) {
if(ncol(cm[[i]])==1){
@@ -23,7 +21,7 @@ qi.vglm <- function (object, simpar, x, x1=NULL, y = NULL) {
all.coef <- NULL
for(i in 1:ndim)
all.coef <- c(all.coef, list(simpar[,v[[i]]]))
- if (model=="mlogit" || model=="oprobit"){
+ if (model=="mlogit"){
if(is.null(colnames(object at y)))
cnames <- ynames <- seq(1,ndim+1,1)
else
@@ -34,14 +32,6 @@ qi.vglm <- function (object, simpar, x, x1=NULL, y = NULL) {
else if(model=="blogit" || model=="bprobit")
cnames <- c("Pr(Y1=0, Y2=0)", "Pr(Y1=0, Y2=1)",
"Pr(Y1=1, Y2=0)", "Pr(Y1=1, Y2=1)")
- else if(model=="sur"){
- if(is.null(colnames(object at y)))
- cnames <- ynames <- seq(1,ndim,1)
- else
- ynames <- cnames <- colnames(object at y)
- for(i in 1:ndim)
- cnames[i] <- paste("E(Y_", ynames[i],"|X)",sep="")
- }
else
stop(paste(model, "is not supported"))
pp.vglm <- function(object, cm, all.coef, x, ndim, cnames){
@@ -59,7 +49,7 @@ qi.vglm <- function (object, simpar, x, x1=NULL, y = NULL) {
}
pr.vglm <- function(object, ev, ynames) { # To assign predicted values.
model <- object$zelig
- if (model == "mlogit" || model == "oprobit") {
+ if (model == "mlogit") {
k <- ncol(ev)
Ipr <- sim.cut <- matrix(NA, nrow = nrow(ev), ncol = ncol(ev))
colnames(Ipr) <- colnames(sim.cut) <- colnames(ev)
@@ -90,9 +80,6 @@ qi.vglm <- function (object, simpar, x, x1=NULL, y = NULL) {
colnames(pr) <- c("(Y1=0, Y2=0)", "(Y1=0, Y2=1)", "(Y1=1, Y2=0)",
"(Y1=1, Y2=1)")
}
-# else if (model == "sur") { # Need to test/add pr.
-# pr[i,] <- mvrnorm(length(ev[i,], mu = ev[i,], vcov = vcov(object)))
-# }
pr
}
if (nrow(x) == 1) {
@@ -100,13 +87,11 @@ qi.vglm <- function (object, simpar, x, x1=NULL, y = NULL) {
pr <- pr.vglm(object, ev, ynames)
}
else {
- ev <- array(dim = c(nrow(simpar), if(model=="sur") ndim else (ndim+1), nrow(x)))
- if (model == "mlogit" || model == "oprobit")
+ ev <- array(dim = c(nrow(simpar), ndim+1, nrow(x)))
+ if (model == "mlogit")
pr <- matrix(nrow=nrow(simpar), ncol=nrow(x))
else if (model == "blogit" || model == "bprobit")
pr <- array(dim = c(nrow(simpar), 4, nrow(x)))
- else if (model == "sur") #### Need to change this
- pr <- NULL
for (i in 1:nrow(x)){
tmp <- matrix(x[i,], nrow=1)
colnames(tmp) <- colnames(x)
@@ -116,20 +101,16 @@ qi.vglm <- function (object, simpar, x, x1=NULL, y = NULL) {
ev[,,i] <- tmp.ev
if (model == "blogit" || model == "bprobit")
pr[,,i] <- tmp.pr
- else if (model == "mlogit" || model == "oprobit")
+ else if (model == "mlogit")
pr[,i] <- tmp.pr
- else if (model == "sur") ##### Placeholder NEED TO FIX!!!!
- pr <- NULL
}
dimnames(ev) <- list(rownames(tmp.ev), colnames(tmp.ev), NULL)
if (model == "blogit" || model == "bprobit")
dimnames(pr) <- list(rownames(tmp.pr), colnames(tmp.pr), NULL)
- else if (model == "mlogit" || model == "oprobit")
+ else if (model == "mlogit")
dimnames(pr) <- list(c(1:nrow(simpar)), NULL)
- else if (model == "sur") ### Need to change.
- dimnames(pr) <- NULL
}
- if (model=="mlogit" || model=="oprobit") {
+ if (model=="mlogit") {
qi <- list(ev=ev, pr=pr)
qi.name <- list(ev="Predicted Probabilities: Pr(Y=k|X)",
pr="Predicted Values: Y=k|X")
@@ -139,16 +120,11 @@ qi.vglm <- function (object, simpar, x, x1=NULL, y = NULL) {
qi.name <- list(ev="Predicted Probabilities: Pr(Y1=k,Y2=l|X)",
pr="Predicted Values: (Y1,Y2)|X")
}
- else if (model=="sur") {
- qi <- list(ev=ev, pr=pr)
- qi.name <- list(ev="Expected Values: E(Yj|X)",
- pr="Predicted Values: Yj|X")
- }
if (!is.null(x1)) {
if (nrow(x1) == 1)
ev1 <- pp.vglm(object, cm, all.coef, x1, ndim, cnames)
else {
- ev1 <- array(dim = c(nrow(simpar), if(model=="sur") ndim else (ndim+1), nrow(x1)))
+ ev1 <- array(dim = c(nrow(simpar), ndim+1, nrow(x1)))
for (i in 1:nrow(x1)) {
tmp <- matrix(x1[i,], nrow=1)
colnames(tmp) <- colnames(x1)
@@ -159,7 +135,7 @@ qi.vglm <- function (object, simpar, x, x1=NULL, y = NULL) {
dimnames(ev1) <- list(rownames(tmp), colnames(tmp), NULL)
}
qi$fd <- ev1 - ev
- if (model=="mlogit" || model=="oprobit") {
+ if (model=="mlogit") {
qi$rr <- ev1 / ev
qi.name$fd <- "First Differences: Pr(Y=k|X1) - Pr(Y=k|X)"
qi.name$rr <- "Risk Ratios: Pr(Y=k|X1) / Pr(Y=k|X)"
@@ -169,8 +145,6 @@ qi.vglm <- function (object, simpar, x, x1=NULL, y = NULL) {
qi.name$fd <- "First Differences: Pr(Y1=k,Y2=l|X1) - Pr(Y1=k,Y2=l|X)"
qi.name$rr <- "Risk Ratios: Pr(Y1=k,Y2=l|X1) / Pr(Y1=k,Y2=l|X)"
}
- else if(model=="sur")
- qi.name$fd <- "First Differences: E(Yj|X1j) - E(Yj|Xj)"
}
if (!is.null(y)) {
tmp.ev <- tmp.pr <- array(NA, dim = dim(qi$ev))
diff --git a/R/relogit.R b/R/relogit.R
index bcb3ef8..24fcb98 100644
--- a/R/relogit.R
+++ b/R/relogit.R
@@ -6,8 +6,8 @@ relogit <- function(formula, data=sys.parent(), tau=NULL,
tau <- unique(tau)
if (length(tau) > 2)
stop("tau must be a vector of length less than or equal to 2")
- mf[[1]] <- as.name("glm")
- mf$family <- as.name("binomial")
+ mf[[1]] <- glm
+ mf$family <- binomial(link="logit")
res <- eval(as.call(mf))
res$call <- match.call(expand.dots = TRUE)
## prior correction
diff --git a/R/setx.MI.R b/R/setx.MI.R
index 4180f2f..7458889 100644
--- a/R/setx.MI.R
+++ b/R/setx.MI.R
@@ -15,7 +15,7 @@ setx.MI <- function(object, fn = list(numeric = mean, ordered =
}
X <- setx(object[[1]], fn = fn, data = dta, cond = FALSE,
counter = NULL, ...)
- class(X) <- c("setx.MI", "setx")
+ class(X) <- c("setx.MI", "setx", "data.frame")
}
else { # conditional prediction
X <- list()
@@ -24,7 +24,7 @@ setx.MI <- function(object, fn = list(numeric = mean, ordered =
for (i in 1:M)
X[[i]] <- setx(object[[i]], fn = NULL, data = data[[i]], cond = TRUE,
counter = counter, ...)
- class(X) <- c("setx.MI", "cond")
+ class(X) <- c("setx.MI", "cond", "data.frame")
}
return(X)
}
diff --git a/R/setx.default.R b/R/setx.default.R
index daef0cc..28694b6 100644
--- a/R/setx.default.R
+++ b/R/setx.default.R
@@ -55,13 +55,13 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
else
tt <- terms(object)
if (is.null(data))
- dta <- eval(object$call$data, sys.parent())
+ odta <- eval(object$call$data, sys.parent())
else
- dta <- data
- data <- dta <- dta[complete.cases(model.frame(tt, dta)), ]
+ odta <- data
+ data <- dta <- model.frame(delete.response(tt), odta)
vars <- names(dta)
if (!is.null(counter)) {
- if (!any(counter == names(model.frame(tt, dta))))
+ if (!any(counter == vars))
stop(paste("the variable specified for counter is not used in the model"))
treat <- data[, names(data)==counter]
if(is.numeric(treat)) {
@@ -95,10 +95,6 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
warning(paste("when cond = TRUE, fn is coerced to NULL"))
fn <- NULL
}
- if (class(object)[1] == "vglm")
- Y <- model.response(model.frame(tt, data=dta), "any")
- else
- Y <- model.extract(model.frame(tt, data=dta), "response")
}
else if (!is.null(fn)) {
if (is.null(fn$numeric) || !is.function(fn$numeric)) {
@@ -120,17 +116,14 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
warning("the only available fn for other is mode.")
fn$other <- mode
}
- for (i in 1:ncol(dta)) {
- v <- na.omit(dta[,i])
- if (any(vars[i] == names(model.frame(tt, dta)))) {
- if (is.numeric(v))
- value <- lapply(list(v), fn$numeric)[[1]]
- else if (is.ordered(v))
- value <- lapply(list(v), fn$ordered)[[1]]
- else
- value <- lapply(list(v), fn$other)[[1]]
- data[,i] <- value
- }
+ for (i in 1:ncol(data)) {
+ if (is.numeric(data[,i]))
+ value <- lapply(list(data[,i]), fn$numeric)[[1]]
+ else if (is.ordered(data[,i]))
+ value <- lapply(list(data[,i]), fn$ordered)[[1]]
+ else
+ value <- lapply(list(data[,i]), fn$other)[[1]]
+ data[,i] <- value
}
opt <- vars[na.omit(pmatch(names(mf), vars))]
maxl <- 1
@@ -146,7 +139,7 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
else
stop("vector inputs should have the same length.")
if (is.factor(data[,opt[i]]))
- data[,opt[i]] <- list(as.factor(value))
+ data[,opt[i]] <- list(factor(value, levels=levels(data[,opt[i]])))
else if (is.numeric(data[,opt[i]]))
data[,opt[i]] <- list(as.numeric(value))
else if (is.logical(data[,opt[i]]))
@@ -156,21 +149,18 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
}
data <- data[1:maxl,]
}
- X <- model.matrix(delete.response(tt), data = data, contrasts =
- if (length(object$contrasts)>0) object$contrasts
- else NULL, xlev = object$xlevels, ...)
- row.names(X) <- 1:nrow(X)
- X <- as.matrix(X)
- class(X) <- "setx"
if (cond) {
- X <- cbind(Y, X)
+ X <- model.frame(tt, odta)
if (!is.null(counter)) {
X <- list(treat=X[treat==1,], control=X[treat==0,])
- class(X$treat) <- class(X$control) <- c("setx", "cond")
+ class(X$treat) <- class(X$control) <- c("data.frame", "cond")
class(X) <- "setx.counter"
}
else
- class(X) <- c("setx", "cond")
+ class(X) <- c("data.frame", "cond")
}
+ else
+ X <- as.data.frame(model.matrix(delete.response(tt), data = data))
+
return(X)
}
diff --git a/R/setx.strata.R b/R/setx.strata.R
index ca02636..bdd8749 100644
--- a/R/setx.strata.R
+++ b/R/setx.strata.R
@@ -38,7 +38,7 @@ setx.strata <- function(object, fn = list(numeric = mean, ordered =
}
}
names(x) <- names(object)
- class(x) <- "setx.strata"
+ class(x) <- c("setx.strata", "data.frame")
return(x)
}
diff --git a/R/sim.cond.R b/R/sim.cond.R
index 376383e..86862cf 100644
--- a/R/sim.cond.R
+++ b/R/sim.cond.R
@@ -5,22 +5,10 @@ sim.cond <- function(object, x, x1=NULL, num=c(1000, 100),
if (!is.null(x1)) {
warning("First Differences are not calculated in conditional prediction models.")
x1 <- NULL
- }
- if (object$call$model %in% c("bprobit", "blogit")) {
- yvar <- x[, 1:2]
- x <- x[, 3:ncol(x)]
}
- else if (any(class(object) == "EI")) {
- yvar <- x[, 1:object$dims[2]]
- x <- x[, (object$dims[2]+1):ncol(x)]
- }
- else if (any(class(object) == "latent"))
- yvar <- NULL
- else {
- yvar <- x[,1]
- x <- x[,2:ncol(x)]
- }
- class(x) <- c("cond")
+ xvar <- model.matrix(object, data = x)
+ yvar <- model.response(x)
+ class(xvar) <- c("matrix", "cond")
if (any(class(object) == "MCMCZelig"))
num <- nrow(object$coefficients)
if (length(num) == 2) {
@@ -33,7 +21,7 @@ sim.cond <- function(object, x, x1=NULL, num=c(1000, 100),
if (!bootstrap & any(class(object) != "relogit"))
simpar <- param(object, num=num, bootstrap=bootstrap)
else if (any(class(object) == "relogit"))
- simpar <- param.relogit(object, num=num, x=x,
+ simpar <- param.relogit(object, num=num, x=xvar,
bootstrap=bootstrap, bootfn=bootfn, ...)
else {
tt <- terms(object)
@@ -53,14 +41,14 @@ sim.cond <- function(object, x, x1=NULL, num=c(1000, 100),
simpar <- prev
}
if (class(object)[1] == "survreg")
- simqi <- qi(object, simpar = simpar, x = x, x1 = x1, y = yvar,
+ simqi <- qi(object, simpar = simpar, x = xvar, x1 = x1, y = yvar,
cond.data = cond.data)
else
- simqi <- qi(object, simpar = simpar, x = x, x1 = x1, y = yvar)
- class(x) <- "cond"
+ 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=x, x1=x1, call = c, zelig.call = object$call,
+ res <- list(x=xvar, x1=x1, call = c, 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 9af8614..876621b 100644
--- a/R/sim.default.R
+++ b/R/sim.default.R
@@ -1,6 +1,9 @@
sim.default <- function(object, x, x1=NULL, num=c(1000, 100),
prev = NULL, bootstrap = FALSE, bootfn=NULL,
cond.data = NULL, ...) {
+ x <- as.matrix(x)
+ if (!is.null(x1))
+ x1 <- as.matrix(x1)
if (any(class(object) == "MCMCZelig"))
num <- nrow(object$coefficients)
else if (length(num) == 2) {
diff --git a/R/summary.glm.robust.R b/R/summary.glm.robust.R
index 5a1ba39..8556a98 100644
--- a/R/summary.glm.robust.R
+++ b/R/summary.glm.robust.R
@@ -1,10 +1,13 @@
summary.glm.robust <- function(object, ...) {
class(object) <- "glm"
res <- summary.glm(object, ...)
- if (is.null(object$robust))
+ if (is.null(object$robust)) {
res$cov.unscaled <- covmat.unscaled <- vcovHAC(object)
+ res$robust <- "vcovHAC"
+ }
else {
fn <- object$robust$method
+ res$robust <- object$robust$method
object$robust$method <- NULL
arg <- object$list
arg$x <- object
@@ -25,5 +28,6 @@ summary.glm.robust <- function(object, ...) {
else
res$coefficients[,4] <- 2 * pt(-abs(tvalue), object$df.residual)
}
+ class(res) <- c("summary.glm.robust","summary.glm")
return(res)
}
diff --git a/R/summary.lm.robust.R b/R/summary.lm.robust.R
index 40126a5..b55b09f 100644
--- a/R/summary.lm.robust.R
+++ b/R/summary.lm.robust.R
@@ -1,12 +1,15 @@
summary.lm.robust <- function(object, ...) {
class(object) <- "lm"
res <- summary.lm(object, ...)
- if (is.null(object$robust))
+ if (is.null(object$robust)) {
res$cov.unscaled <- R <- vcovHC(object)/(res$sigma^2)
+ res$robust <- "vcovHC"
+ }
else {
fn <- object$robust$method
+ res$robust <- object$robust$method
object$robust$method <- NULL
- arg <- object$list
+ arg <- object$robust
arg$x <- object
res$cov.unscaled <- R <- eval(do.call(fn, arg))/(res$sigma^2)
}
@@ -16,6 +19,8 @@ summary.lm.robust <- function(object, ...) {
dimnames(res$correlation) <- dimnames(res$cov.unscaled)
}
res$coefficients[,3] <- tval <- coefficients(object)/se
- res$coefficients[,4] <- 2*pt(abs(tval), res$df[2], lower.tail = FALSE)
+ res$coefficients[,4] <- 2*pt(abs(tval), res$df[2], lower.tail =
+ FALSE)
+ class(res) <- c("summary.lm.robust", "summary.lm")
return(res)
}
diff --git a/R/zelig.R b/R/zelig.R
index 3f211c2..d8b90f6 100644
--- a/R/zelig.R
+++ b/R/zelig.R
@@ -1,9 +1,9 @@
-zelig <- function(formula, model, data, by = NULL, robust = FALSE, ...) {
- fn <- paste("zelig2", model, sep = "")
- if (!exists(fn))
+zelig <- function(formula, model, data, by = NULL, ...) {
+ fn1 <- paste("zelig2", model, sep = "")
+ fn2 <- paste("zelig3", model, sep = "")
+ if (!exists(fn1))
stop(model, "not supported. Type help.zelig(\"models\") to list supported models.")
- mf <- match.call(expand.dots = TRUE)
- mf$robust <- NULL
+ mf <- zelig.call <- match.call(expand.dots = TRUE)
if (missing(by))
by <- NULL
N <- M <- 1
@@ -22,7 +22,7 @@ zelig <- function(formula, model, data, by = NULL, robust = FALSE, ...) {
lev <- sort(unique(idx))
N <- length(lev)
}
- mf <- do.call(fn, list(formula, model, dat, N, ...))
+ mf <- do.call(fn1, list(formula, model, dat, N, ...))
for (i in 1:N) {
if (N > 1) {
dat <- list()
@@ -45,32 +45,12 @@ zelig <- function(formula, model, data, by = NULL, robust = FALSE, ...) {
d <- d[complete.cases(model.frame(as.formula(formula), d)),]
mf$data <- d
res <- eval(as.call(mf))
- res$call <- match.call()
+ if (exists(fn2))
+ res <- do.call(fn2, list(res = res, fcall = as.list(mf),
+ zcall = as.list(zelig.call)))
+ res$call <- zelig.call
res$data <- res$call$data
res$zelig <- model
- if (is.list(robust)) {
- if (any(c("lm", "glm") %in% class(res)[1])) {
- require(sandwich)
- ctmp <- class(res)
- class(res) <- c(paste(ctmp[1], ".robust", sep=""), ctmp)
- if (!any(robust$method %in% c("vcovHC", "vcovHAC", "kernHAC")))
- stop("such a robust option is not supported")
- else if ((robust$method == "vcovHC") & ("lm" != class(res)[1]))
- stop("vcovHC is supported only for ols")
- res$robust <- robust
- }
- else
- stop("robust option is not supported for this model.")
- }
- else if (robust) {
- if (any(c("lm", "glm") %in% class(res)[1])) {
- require(sandwich)
- ctmp <- class(res)
- class(res) <- c(paste(ctmp[1], ".robust", sep=""), ctmp)
- }
- else
- stop("robust option is not supported for this model.")
- }
if (M > 1)
obj[[j]] <- res
else
diff --git a/R/zelig2beta.R b/R/zelig2beta.R
index cbd5836..4e4d91e 100644
--- a/R/zelig2beta.R
+++ b/R/zelig2beta.R
@@ -1,6 +1,6 @@
zelig2beta <- 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]] <- as.name("beta.regression")
as.call(mf)
}
diff --git a/R/zelig2blogit.R b/R/zelig2blogit.R
index 2e42c4a..63fe184 100644
--- a/R/zelig2blogit.R
+++ b/R/zelig2blogit.R
@@ -6,7 +6,7 @@ zelig2blogit <- function(formula, model, data, M, constrain = NULL,
else
stop("Please install VGAM using \n install.packages(\"VGAM\", CRAN = \"http://www.stat.auckland.ac.nz/~yee\")")
mf <- match.call(expand.dots = TRUE)
- mf[[1]] <- as.name("vglm")
+ mf[[1]] <- VGAM::vglm
mf$family <- as.name("blogit")
mf$... <- NULL
tmp <- cmvglm(formula, model, constrain, omit, constant, 3)
@@ -14,5 +14,6 @@ zelig2blogit <- function(formula, model, data, M, constrain = NULL,
mf$constraints <- tmp$constraints
blogit <<- function() binom2.or(zero=NULL)
mf$model <- mf$constrain <- mf$omit <- mf$constant <- mf$M <- NULL
+ mf$robust <- NULL
as.call(mf)
}
diff --git a/R/zelig2bprobit.R b/R/zelig2bprobit.R
index 25baf89..2c28116 100644
--- a/R/zelig2bprobit.R
+++ b/R/zelig2bprobit.R
@@ -6,7 +6,7 @@ zelig2bprobit <- function(formula, model, data, M, constrain = NULL,
else
stop("Please install VGAM using \n install.packages(\"VGAM\", CRAN = \"http://www.stat.auckland.ac.nz/~yee\")")
mf <- match.call(expand.dots = TRUE)
- mf[[1]] <- as.name("vglm")
+ mf[[1]] <- VGAM::vglm
mf$family <- as.name("bprobit")
tmp <- cmvglm(formula, model, constrain, omit, constant, 3)
mf$formula <- tmp$formula
diff --git a/R/zelig2exp.R b/R/zelig2exp.R
index 7a36c62..d8becda 100644
--- a/R/zelig2exp.R
+++ b/R/zelig2exp.R
@@ -2,7 +2,24 @@ zelig2exp <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
mf$model <- mf$M <- NULL
require(survival)
- mf[[1]] <- as.name("survreg")
- mf$dist <- as.character("exponential")
+ mf[[1]] <- survival::survreg
+ mf$dist <- "exponential"
+ 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("+", " cluster(",
+ mf$cluster, ")")))
+ mf$cluster <- NULL
+ }
+ else if (mf$robust)
+ mf$formula <- as.formula(paste(paste(deparse(formula[[2]])),
+ paste(deparse(formula[[1]])),
+ paste(deparse(formula[[3]])),
+ paste("+", " cluster(1:nrow(",
+ deparse(formula[[2]]),
+ "))")))
as.call(mf)
}
diff --git a/R/zelig2gamma.R b/R/zelig2gamma.R
index f9150a7..afd3ff6 100644
--- a/R/zelig2gamma.R
+++ b/R/zelig2gamma.R
@@ -1,7 +1,7 @@
zelig2gamma <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- NULL
- mf[[1]] <- as.name("glm")
- mf$family <- as.name("Gamma")
+ mf$model <- mf$M <- mf$robust <- NULL
+ mf[[1]] <- glm
+ mf$family <- Gamma
as.call(mf)
}
diff --git a/R/zelig2logit.R b/R/zelig2logit.R
index 7cd2612..d176d79 100644
--- a/R/zelig2logit.R
+++ b/R/zelig2logit.R
@@ -1,7 +1,7 @@
zelig2logit <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- NULL
- mf[[1]] <- as.name("glm")
- mf$family <- as.name("binomial")
+ mf$model <- mf$M <- mf$robust <- NULL
+ mf[[1]] <- glm
+ mf$family <- binomial(link="logit")
as.call(mf)
}
diff --git a/R/zelig2lognorm.R b/R/zelig2lognorm.R
index b8f3714..2c56a0f 100644
--- a/R/zelig2lognorm.R
+++ b/R/zelig2lognorm.R
@@ -2,7 +2,24 @@ zelig2lognorm <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
mf$model <- mf$M <- NULL
require(survival)
- mf[[1]] <- as.name("survreg")
- mf$dist <- as.character("lognormal")
+ mf[[1]] <- survival::survreg
+ mf$dist <- "lognormal"
+ 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("+", " cluster(",
+ mf$cluster, ")")))
+ mf$cluster <- NULL
+ }
+ else if (mf$robust)
+ mf$formula <- as.formula(paste(paste(deparse(formula[[2]])),
+ paste(deparse(formula[[1]])),
+ paste(deparse(formula[[3]])),
+ paste("+", " cluster(1:nrow(",
+ deparse(formula[[2]]),
+ "))")))
as.call(mf)
}
diff --git a/R/zelig2ls.R b/R/zelig2ls.R
index 8e4b22f..e2af63f 100644
--- a/R/zelig2ls.R
+++ b/R/zelig2ls.R
@@ -1,6 +1,6 @@
zelig2ls <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- NULL
- mf[[1]] <- as.name("lm")
+ mf$model <- mf$M <- mf$robust <- NULL
+ mf[[1]] <- lm
as.call(mf)
}
diff --git a/R/zelig2mlogit.R b/R/zelig2mlogit.R
index 0f79ce1..871e811 100644
--- a/R/zelig2mlogit.R
+++ b/R/zelig2mlogit.R
@@ -5,8 +5,8 @@ zelig2mlogit <- function(formula, model, data, M, ...) {
else
stop("Please install VGAM using \n install.packages(\"VGAM\", CRAN = \"http://www.stat.auckland.ac.nz/~yee\")")
mf <- match.call(expand.dots = TRUE)
- mf[[1]] <- as.name("vglm")
- mf$family <- as.name("multinomial")
+ mf[[1]] <- VGAM::vglm
+ mf$family <- VGAM::multinomial
if (mf$M == 1)
ndim <- length(unique(na.omit((eval(mf$formula[[2]], mf$data))))) - 1
if (mf$M > 1) {
diff --git a/R/zelig2mloglm.R b/R/zelig2mloglm.R
index 3068456..94141cc 100644
--- a/R/zelig2mloglm.R
+++ b/R/zelig2mloglm.R
@@ -1,7 +1,7 @@
zelig2mloglm <- function(formula, model, data, M, ...) {
mf$model <- mf$M <- mf$... <- NULL
require(nnet)
- mf[[1]] <- as.name("multinom")
+ mf[[1]] <- nnet::multinom
mf$Hess <- TRUE
as.call(mf)
}
diff --git a/R/zelig2negbin.R b/R/zelig2negbin.R
index 6725da6..c4cff1d 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[[1]] <- as.name("glm.nb")
+ mf[[1]] <- MASS::glm.nb
as.call(mf)
}
diff --git a/R/zelig2normal.R b/R/zelig2normal.R
index f1a3bb5..5a9f884 100644
--- a/R/zelig2normal.R
+++ b/R/zelig2normal.R
@@ -1,7 +1,7 @@
zelig2normal <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- NULL
- mf[[1]] <- as.name("glm")
- mf$family <- as.name("gaussian")
+ mf$model <- mf$M <- mf$robust <- NULL
+ mf[[1]] <- glm
+ mf$family <- gaussian
as.call(mf)
}
diff --git a/R/zelig2ologit.R b/R/zelig2ologit.R
index 2360654..82509d3 100644
--- a/R/zelig2ologit.R
+++ b/R/zelig2ologit.R
@@ -1,7 +1,8 @@
zelig2ologit <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
mf$model <- mf$M <- NULL
- mf[[1]] <- as.name("polr")
+ mf[[1]] <- MASS::polr
mf$Hess <- TRUE
+ mf$method <- "logistic"
as.call(mf)
}
diff --git a/R/zelig2oprobit.R b/R/zelig2oprobit.R
index d50a89e..410f472 100644
--- a/R/zelig2oprobit.R
+++ b/R/zelig2oprobit.R
@@ -1,27 +1,8 @@
zelig2oprobit <- function(formula, model, data, M, ...) {
- check <- library()
- if(any(check$results[,"Package"] == "VGAM"))
- require(VGAM)
- else
- stop("Please install VGAM using \n install.packages(\"VGAM\", CRAN = \"http://www.stat.auckland.ac.nz/~yee\")")
mf <- match.call(expand.dots = TRUE)
- if(is.null(mf$equal))
- ordinal.probit <<- function() cumulative(link=probit, parallel=T)
- else
- ordinal.probit <<- function() cumulative(link=probit)
- mf[[1]] <- as.name("vglm")
- mf$family <- as.name("ordinal.probit")
- if (mf$M == 1)
- ndim <- length(unique(na.omit((eval(mf$formula[[2]], mf$data))))) - 1
- if (mf$M > 1) {
- ntmp <- array()
- for (i in 1:mf$M)
- ntmp[i] <- length(unique(na.omit((eval(mf$formula[[2]], mf$data[[i]]))))) - 1
- ndim <- max(ntmp)
- }
- tmp <- cmvglm(mf$formula, mf$model, mf$equal, mf$zeros, mf$ones, ndim)
- mf$formula <- tmp$formula
- mf$constraints <- tmp$constraints
- mf$model <- mf$equal <- mf$ones <- mf$zeros <- mf$M <- NULL
+ mf$model <- mf$M <- NULL
+ mf[[1]] <- MASS::polr
+ mf$Hess <- TRUE
+ mf$method <- "probit"
as.call(mf)
}
diff --git a/R/zelig2poisson.R b/R/zelig2poisson.R
index 4edabba..11b4e2e 100644
--- a/R/zelig2poisson.R
+++ b/R/zelig2poisson.R
@@ -1,7 +1,7 @@
zelig2poisson <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- NULL
- mf[[1]] <- as.name("glm")
- mf$family <- as.name("poisson")
+ mf$model <- mf$M <- mf$robust <- NULL
+ mf[[1]] <- glm
+ mf$family <- poisson
as.call(mf)
}
diff --git a/R/zelig2probit.R b/R/zelig2probit.R
index be5a03e..9c35621 100644
--- a/R/zelig2probit.R
+++ b/R/zelig2probit.R
@@ -1,8 +1,7 @@
zelig2probit <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- NULL
- mf[[1]] <- as.name("glm")
- mf$family <- as.name("binomial.probit")
- binomial.probit <<- function() binomial(link="probit")
+ mf$model <- mf$M <- mf$probit <- NULL
+ mf[[1]] <- stats::glm
+ mf$family <- binomial(link="probit")
as.call(mf)
}
diff --git a/R/zelig2weibull.R b/R/zelig2weibull.R
index a633439..75e46e1 100644
--- a/R/zelig2weibull.R
+++ b/R/zelig2weibull.R
@@ -2,7 +2,24 @@ zelig2weibull <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
mf$model <- mf$M <- NULL
require(survival)
- mf[[1]] <- as.name("survreg")
- mf$dist <- as.character("weibull")
+ mf[[1]] <- survival::survreg
+ mf$dist <- "weibull"
+ 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("+", " cluster(",
+ mf$cluster, ")")))
+ mf$cluster <- NULL
+ }
+ else if (mf$robust)
+ mf$formula <- as.formula(paste(paste(deparse(formula[[2]])),
+ paste(deparse(formula[[1]])),
+ paste(deparse(formula[[3]])),
+ paste("+", " cluster(1:nrow(",
+ deparse(formula[[2]]),
+ "))")))
as.call(mf)
}
diff --git a/R/zelig3glm.R b/R/zelig3glm.R
new file mode 100644
index 0000000..3349769
--- /dev/null
+++ b/R/zelig3glm.R
@@ -0,0 +1,22 @@
+zelig3logit <- zelig3probit <- zelig3normal <- zelig3gamma <-
+ zelig3poisson <- zelig3negbin <- zelig3relogit <-
+ function(res, fcall = NULL, zcall = NULL) {
+ rob <- zcall$robust
+ if (!is.null(rob)) {
+ require(sandwich)
+ if (is.list(rob)) {
+ if (!any(rob$method %in% c("vcovHAC", "kernHAC", "weave")))
+ stop("such a robust option is not supported")
+ else {
+ class(res) <- c("glm.robust", class(res))
+ res$robust <- rob
+ }
+ }
+ else if (!is.logical(rob))
+ stop("invalid input for robust. Choose either TRUE or a list of options.")
+ else if (rob)
+ class(res) <- c("glm.robust", class(res))
+ }
+ return(res)
+ }
+
diff --git a/R/zelig3ls.R b/R/zelig3ls.R
new file mode 100644
index 0000000..db5662b
--- /dev/null
+++ b/R/zelig3ls.R
@@ -0,0 +1,19 @@
+zelig3ls <- function(res, fcall = NULL, zcall = NULL) {
+ rob <- eval(zcall$robust)
+ if (!is.null(rob)) {
+ require(sandwich)
+ if (is.list(rob)) {
+ if (!any(rob$method %in% c("vcovHC", "vcovHAC",
+ "kernHAC", "weave")))
+ stop("such a robust option is not supported")
+ res$robust <- rob
+ class(res) <- c("lm.robust", class(res))
+ }
+ else if (!is.logical(rob))
+ stop("invalid input for robust. Choose either TRUE or a list of options.")
+ else if (rob)
+ class(res) <- c("lm.robust", class(res))
+ }
+ return(res)
+}
+
diff --git a/R/zelig3ologit.R b/R/zelig3ologit.R
new file mode 100644
index 0000000..f835c06
--- /dev/null
+++ b/R/zelig3ologit.R
@@ -0,0 +1,12 @@
+zelig3ologit <- function(res, fcall = NULL, zcall = NULL) {
+ inv.link <- function(eta, zeta) {
+ tmp1 <- matrix(1, nrow = length(eta), ncol = length(zeta) + 1)
+ ilogit <- function(e, z) {
+ exp(z - e) / (1 + exp(z - e))
+ }
+ tmp1[, 1:length(zeta)] <- sapply(zeta, ilogit, e = eta)
+ tmp1
+ }
+ res$inv.link <- as.function(inv.link)
+ res
+}
diff --git a/R/zelig3oprobit.R b/R/zelig3oprobit.R
new file mode 100644
index 0000000..157823b
--- /dev/null
+++ b/R/zelig3oprobit.R
@@ -0,0 +1,12 @@
+zelig3oprobit <- function(res, fcall = NULL, zcall = NULL) {
+ inv.link <- function(eta, zeta) {
+ tmp1 <- matrix(1, nrow = length(eta), ncol = length(zeta) + 1)
+ iprobit <- function(z, e)
+ pnorm(z - e)
+ tmp1[, 1:length(zeta)] <- sapply(zeta, iprobit, e = eta)
+ tmp1
+ }
+ res$inv.link <- as.function(inv.link)
+ res
+}
+
diff --git a/README b/README
index 985803e..6a9ef87 100644
--- a/README
+++ b/README
@@ -1,6 +1,12 @@
+2.2-1 (July 11, 2005): Stable release for R 2.0.0-2.1.0. Revised
+ ordinal probit to use MASS library. Added robust standard errors
+ for the following regression models: exponential, gamma, logit,
+ lognormal, least squares, negative binomial, normal (Gaussian),
+ poisson, probit, and weibull.
+
2.1-4 (May 22, 2005): Stable release for R 1.9.1-2.1.0. Revised
- help.zelig() to deal with CRAN build of Windows version. Added
- recode of slots to lists in NAMESPACE. Revised install.R script
+ help.zelig() to deal with CRAN build of Windows version. Added
+ recode of slots to lists in NAMESPACE. Revised install.R script
to deal with changes to install.packages(). (reported by Dan
Powers and Ying Lu)
diff --git a/demo/00Index b/demo/00Index
index 7ee10ce..280023d 100644
--- a/demo/00Index
+++ b/demo/00Index
@@ -17,6 +17,7 @@ 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
diff --git a/demo/oprobit.R b/demo/oprobit.R
index 4dad7f1..b5bfc82 100644
--- a/demo/oprobit.R
+++ b/demo/oprobit.R
@@ -8,7 +8,7 @@ user.prompt()
z.out1 <- zelig(as.factor(cost) ~ mil + coop, model = "oprobit",
data = sanction)
user.prompt()
-print(summary(z.out1))
+summary(z.out1)
# Set the explanatory variables to their means, with 'mil' set
# to 0 (no military action in addition to sanctions) in the baseline
@@ -23,8 +23,7 @@ x.high <- setx(z.out1, mil = 1)
user.prompt()
s.out1 <- sim(z.out1, x = x.low, x1 = x.high)
user.prompt()
-print(summary(s.out1))
-
+summary(s.out1)
##### Example 2: Creating An Ordered Dependent Variable #####
@@ -38,7 +37,7 @@ sanction$ncost <- factor(sanction$ncost, ordered = TRUE,
user.prompt()
z.out2 <- zelig(ncost ~ mil + coop, model = "oprobit", data = sanction)
user.prompt()
-print(summary(z.out2))
+summary(z.out2)
# Set the explanatory variables to their observed values:
user.prompt()
@@ -48,7 +47,7 @@ x.out2 <- setx(z.out2, fn = NULL)
user.prompt()
s.out2 <- sim(z.out2, x = x.out2)
user.prompt()
-print(summary(s.out2))
+summary(s.out2)
diff --git a/demo/robust.R b/demo/robust.R
index 78d8288..a816aed 100644
--- a/demo/robust.R
+++ b/demo/robust.R
@@ -30,8 +30,37 @@ plot(s.out1)
data(hoff)
# Fit the model with robust standard error
user.prompt()
-z.out1 <- zelig(L2SocSec ~ Just503D + Just503R + Just503D:RGovDumy +
+z.out2 <- zelig(L2SocSec ~ Just503D + Just503R + Just503D:RGovDumy +
Just503R:I(1-RGovDumy), model = "ls", data = hoff,
robust = list(method="vcovHAC", order.by=hoff$year, adjust=TRUE))
user.prompt()
-print(summary(z.out1))
+print(summary(z.out2))
+
+##### Example 3: weibull regression with
+##### heteroskedasticity consistent standard errors
+##### and using invest as a cluster
+
+# Attach sample data and variable names:
+data(coalition)
+# Fit the model with robust standard error
+user.prompt()
+z.out3 <- zelig(Surv(duration, ciep12) ~ polar + numst2 +
+ crisis, model = "weibull", data = coalition,
+ cluster = "invest", robust = TRUE)
+user.prompt()
+print(summary(z.out3))
+
+
+#####
+##### Example 4: logit regression with heteroskedasticity and
+##### autocorrelation consistent standard errors
+
+# Attach sample data and variable names
+data(turnout)
+# Fit the model with robust standrad error
+user.prompt()
+z.out4 <- zelig(vote ~ race + educate, model = "logit",
+ data = turnout, robust=TRUE)
+user.prompt()
+print(summary(z.out4))
+
diff --git a/demo/vertci.R b/demo/vertci.R
index 144eb1b..cb5033c 100644
--- a/demo/vertci.R
+++ b/demo/vertci.R
@@ -12,10 +12,8 @@ summary(z.out)
## Creating setx structures with education set to high school and
## post-college levels, for the whole range of the age variable.
user.prompt()
-age.range <- 18:95
-user.prompt()
-x.low <- setx(z.out, educate = 12, age = age.range)
-x.high <- setx(z.out, educate = 16, age = age.range)
+x.low <- setx(z.out, educate = 12, age = 18:95)
+x.high <- setx(z.out, educate = 16, age = 18:95)
## Using sim to generate the simulated predicted probabilites:
user.prompt()
diff --git a/man/help.zelig.Rd b/man/help.zelig.Rd
index 9b5c8d4..4af2258 100644
--- a/man/help.zelig.Rd
+++ b/man/help.zelig.Rd
@@ -11,13 +11,13 @@
}
\usage{
-help.zelig(object)
+help.zelig(...)
}
\arguments{
- \item{object}{a character string representing a Zelig command or model.
- \code{help.zelig("command")} will take you to an index of Zelig
- commands and \code{help.zelig("model")} will take you to a list of
+ \item{...}{a Zelig command or model.
+ \code{help.zelig(command)} will take you to an index of Zelig
+ commands and \code{help.zelig(model)} will take you to a list of
models. }
}
@@ -25,12 +25,6 @@ help.zelig(object)
\url{http://gking.harvard.edu/zelig}.
}
-\references{Gary King, Michael Tomz, and Jason Wittenberg (2000),
-``Making the Most of Statistical Analyses: Improving Interpretation and
- Presentation,'' \emph{American Journal of Political Science}, vol. 44,
- pp.347--355.
-}
-
\author{
Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
<\email{king at harvard.edu}>; Olivia Lau<\email{olau at fas.harvard.edu}>
diff --git a/man/summary.zelig.Rd b/man/summary.zelig.Rd
index 91f41ba..532b6f2 100644
--- a/man/summary.zelig.Rd
+++ b/man/summary.zelig.Rd
@@ -10,7 +10,7 @@
interst.}
\usage{
-\method{summary}{zelig}(object, subset = NULL, CI = 95, stats = c("mean", "sd", "min", "max"), \dots)
+\method{summary}{zelig}(object, subset = NULL, CI = 95, stats = c("mean", "sd"), \dots)
}
\arguments{
--
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