Thanks! I'll get this into including this into the next release.
On Mar 2, 2012, at 3:16 PM, Rogério Barbosa wrote:
Hey all again,
I would like to inform that the most updated version of plot.ci function doesn't
have the item "logit.survey" listed on its algorithm. So it is not possible to
plot confidence intervals for predicted values calculated using this model.
I tried to make that inclusion myself -- and it worked.
Here it goes the modified function.
(I found the original at:
http://r.iq.harvard.edu/src/contrib/Zelig/R/plot.ci.R)
Best,
plot.ci <- function(x, CI=95, qi = "ev", main = "",
ylab = NULL, xlab = NULL, xlim = NULL,
ylim = NULL, col = c("red", "blue"), ...) {
"%w/o%" <- function(x,y) x[!x %in% y] #-- x without y
if (class(x) != "zelig")
stop(" plot.ci() works only for sim() output.")
if (!(x$zelig.call$model) %in%
c("ls", "logit", "logit.survey",
"normal.survey","probit", "exp", "gamma",
"lognorm",
"weibull", "normal", "poisson",
"tobit", "relogit",
"negbin", "logit.bayes", "probit.bayes",
"poisson.bayes", "normal.bayes", "tobit.bayes",
"ls.mixed", "logit.mixed", "probit.mixed",
"gamma.mixed", "poisson.mixed",
"logit.gam", "gamma.gee", "normal.gam",
"poisson.gam",
"probit.gam", "logit.gee", "normal.gee",
"poisson.gee", "probit.gee", "normal.survey"))
stop("\n plot.ci() is valid only for non-categorical, univariate response
models.")
cip <- c((100-CI)/200, 1-(100-CI)/200)
summarize <- function(z, cip){
res <- NULL
res <- cbind(res, apply(z, 2, quantile, prob=cip[1]))
res <- cbind(res, apply(z, 2, quantile, prob=cip[2]))
res
}
vv <- apply(x$x, 2, unique)
idx <- sapply(vv, length)
cidx <- which(idx > 1)
if (!is.null(x$x1)) {
vv1 <- apply(x$x1, 2, unique)
idx1 <- sapply(vv1, length)
cidx1 <- which(idx1 > 1)
if (!identical(names(idx), names(idx1)))
stop("variables in x and x1 do not match.")
## Checking for one dimension of variation, including interaction terms
if (length(cidx) > length(cidx1)) {
tmp <- names(idx)[cidx %w/o% cidx1]
tmp1 <- names(idx)[cidx[cidx %in% cidx1]]
}
else {
tmp <- names(idx1)[cidx1 %w/o% cidx]
tmp1 <- names(idx1)[cidx1[cidx1 %in% cidx]]
}
check <- grep(tmp1, tmp)
if (length(check) != length(tmp))
stop("x and x1 vary on more than one dimension.")
}
var <- vv[[cidx[1]]]
q <- pmatch(qi, names(x$qi))
qofi <- x$qi[[q]]
sum.qi <- summarize(qofi, cip)
if (!is.null(x$x1) && qi == "ev") {
fd <- x$qi$fd
ev1 <- fd + qofi
sum.qi1 <- summarize(ev1, cip)
}
else sum.qi1 <- NULL
if (is.null(ylab)) ylab <- x$qi.name[[q]]
if (is.null(xlab)) xlab <- paste("Range of", colnames(x$x)[cidx[1]])
if (is.null(ylim)) {
if (is.null(sum.qi1)) ylim <- c(min(sum.qi), max(sum.qi))
else ylim <- c(min(sum.qi, sum.qi1), max(sum.qi, sum.qi1))
}
if (is.null(xlim)) xlim <- c(min(var), max(var))
plot.default(var, type = "n", ylab = ylab, main = main, xlab = xlab,
xlim = xlim, ylim = ylim)
for (i in 1:length(var)) {
lines(c(var[i], var[i]), c(sum.qi[i,1], sum.qi[i,2]), col = col[1], ...)
if (!is.null(x$x1) && qi == "ev")
lines(c(var[i], var[i]), c(sum.qi1[i,1], sum.qi1[i,2]), col = col[2], ...)
}
}
- - - - -
Rogério J. Barbosa
Researcher at Centre for Metropolitan Studies/Cebrap
São Paulo - Brazil
-
--
Zelig Mailing List, served by HUIT
Send messages: zelig(a)lists.gking.harvard.edu
[un]subscribe Options:
http://lists.gking.harvard.edu/mailman/listinfo/zelig
Zelig program information:
http://gking.harvard.edu/zelig/
Zelig mailing list
Zelig(a)lists.gking.harvard.edu
https://lists.gking.harvard.edu/mailman/listinfo/zelig