[r-cran-zelig] 18/102: Import Upstream version 2.5-2
Andreas Tille
tille at debian.org
Sun Jan 8 16:58:10 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 9b3a30a2d1a9ba6cec1ae2b9b62c53190cb218ff
Author: Andreas Tille <tille at debian.org>
Date: Sun Jan 8 09:39:02 2017 +0100
Import Upstream version 2.5-2
---
.Rbuildindex.18782 | 0
DESCRIPTION | 8 +-
NAMESPACE | 27 ++-
R/MCMC/class.ind.R | 11 -
R/MCMC/param.MCMCZelig.R | 14 --
R/MCMC/print.summary.MCMCZelig.R | 12 -
R/MCMC/qi.MCMCZelig.R | 393 ---------------------------------
R/MCMC/summary.MCMCZelig.R | 19 --
R/MCMC/zelig2MCMC.R | 348 -----------------------------
R/MCMC/zelig3MCMC.R | 135 -----------
R/MULTIPLE/cmsystemfit.R | 19 --
R/MULTIPLE/model.matrix.multiple.R | 67 ------
R/MULTIPLE/terms.multiple.R | 19 --
R/MULTIPLE/zelig2w2sls.R | 15 --
R/MULTIPLE_NEW/model.frame.list.R | 30 ---
R/MULTIPLE_NEW/model.matrix.list.R | 171 --------------
R/MULTIPLE_NEW/parse.formula.R | 90 --------
R/MULTIPLE_NEW/parse.par.R | 74 -------
R/MULTIPLE_NEW/tag.R | 5 -
R/MULTIPLE_NEW/terms.list.R | 113 ----------
R/MULTIPLE_NEW/toBuildFormula.R | 17 --
R/categories.R | 9 +
R/check.describe.R | 50 +++++
R/cmvglm.R | 71 ++++--
R/current.packages.R | 51 +++++
R/describe.blogit.R | 19 ++
R/describe.bprobit.R | 19 ++
R/describe.default.R | 12 +
R/describe.exp.R | 17 ++
R/describe.factor.bayes.R | 15 ++
R/describe.factor.mix.R | 17 ++
R/describe.factor.ord.R | 17 ++
R/describe.gamma.R | 15 ++
R/describe.irt1d.R | 17 ++
R/describe.irtkd.R | 17 ++
R/describe.logit.R | 15 ++
R/describe.logit.bayes.R | 15 ++
R/describe.lognorm.R | 17 ++
R/describe.ls.R | 15 ++
R/describe.mlogit.R | 18 ++
R/describe.mlogit.bayes.R | 17 ++
R/describe.negbin.R | 15 ++
R/describe.normal.R | 15 ++
R/describe.normal.bayes.R | 15 ++
R/describe.ologit.R | 17 ++
R/describe.oprobit.R | 19 ++
R/describe.oprobit.bayes.R | 17 ++
R/describe.poisson.R | 15 ++
R/describe.poisson.bayes.R | 15 ++
R/describe.probit.R | 15 ++
R/describe.probit.bayes.R | 15 ++
R/describe.relogit.R | 15 ++
R/describe.tobit.bayes.R | 15 ++
R/describe.weibull.R | 17 ++
R/{MULTIPLE_NEW => }/make.parameters.R | 26 ++-
R/model.end.R | 10 +-
R/model.frame.multiple.R | 104 +++++++++
R/model.matrix.multiple.R | 100 +++++++++
R/multi.R | 26 +++
R/parse.formula.R | 323 +++++++++++++++++++++++++++
R/parse.par.R | 50 +++++
R/{MULTIPLE_NEW => }/put.start.R | 3 +-
R/qi.MCMCZelig.R | 18 +-
R/qi.survreg.R | 16 +-
R/relogit.R | 6 +
R/{MULTIPLE_NEW => }/set.start.R | 4 +-
R/setx.default.R | 67 +++---
R/sim.default.R | 1 +
R/sim.setx.MI.R | 3 +-
R/summary.glm.robust.R | 9 +-
R/systemfit/2sls.tex | 1 +
R/systemfit/3sls.tex | 1 +
R/systemfit/sur.tex | 1 +
R/terms.multiple.R | 98 ++++++++
R/ternaryplot.R | 89 ++++++++
R/vdc.R | 244 ++++++++++++++++++++
R/zelig.R | 74 +------
R/{zelig.R => zelig.default.R} | 12 +-
R/zelig2blogit.R | 8 +-
R/zelig2bprobit.R | 11 +-
R/zelig2gamma.R | 3 +-
R/zelig2logit.R | 5 +-
R/zelig2ls.R | 3 +-
R/zelig2mlogit.R | 16 +-
R/zelig2normal.R | 3 +-
R/zelig2poisson.R | 3 +-
R/zelig2probit.R | 5 +-
R/zelig2tobit.R | 59 +++++
R/zelig3glm.R | 2 +-
R/zelig3relogit.R | 40 ++--
README | 21 ++
data/Zelig.url.tab | 2 +
demo/00Index | 79 +++----
demo/bivariate.probit.R | 81 +++++++
demo/blogit.R | 3 +-
demo/bprobit.R | 12 +-
demo/factor.mix.R | 3 +-
demo/factor.ord.R | 4 +-
demo/irt1d.R | 2 +-
demo/match.R | 1 -
demo/mlogit.R | 40 +++-
demo/normal.regression.R | 102 +++++++++
demo/oprobit.bayes.R | 2 +-
demo/tobit.R | 47 ++++
man/current.packages.Rd | 36 +++
man/model.end.Rd | 41 ++++
man/model.frame.multiple.Rd | 60 +++++
man/model.matrix.multiple.Rd | 91 ++++++++
man/parse.formula.Rd | 88 ++++++++
man/parse.par.Rd | 82 +++++++
man/plot.ci.Rd | 2 +-
man/put.start.Rd | 34 +++
man/repl.Rd | 4 +-
man/set.start.Rd | 40 ++++
man/setx.Rd | 12 +-
man/ternaryplot.Rd | 72 ++++++
man/ternarypoints.Rd | 6 +-
man/zeligVDC.Rd | 82 +++++++
118 files changed, 2884 insertions(+), 1844 deletions(-)
diff --git a/.Rbuildindex.18782 b/.Rbuildindex.18782
new file mode 100644
index 0000000..e69de29
diff --git a/DESCRIPTION b/DESCRIPTION
index 59b4dd8..d6ab540 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,13 +1,13 @@
Package: Zelig
-Version: 2.4-5
-Date: 2005-10-18
+Version: 2.5-2
+Date: 2006-02-03
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 (>= 2.0.0), MASS, boot
-Suggests: VGAM, survival, sandwich, zoo, MCMCpack, coda
+Suggests: VGAM, mvtnorm, survival, sandwich (>= 1.1-0), zoo (>= 1.0-4), 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: Tue Oct 18 00:49:19 2005; olau
+Packaged: Sat Feb 4 17:52:31 2006; olau
diff --git a/NAMESPACE b/NAMESPACE
index 7b21a15..f24c91b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -11,17 +11,32 @@ importFrom(stats, quasi)
importFrom(stats, quasibinomial)
importFrom(stats, quasipoisson)
-export( dims,
+export( current.packages,
+ dims,
gsource,
help.zelig,
+ model.end,
+ model.matrix.multiple,
+ model.frame.multiple,
+ parse.formula,
+ parse.par,
plot.ci,
plot.zelig,
+ put.start,
+ repl,
rocplot,
+ set.start,
setx,
sim,
+ ternaryplot,
ternarypoints,
- user.prompt,
- zelig
+ user.prompt,
+ zelig,
+ zeligListModels,
+ zeligInstalledModels,
+ zeligModelDependency,
+ zeligGetSpecial,
+ zeligDescribeModelXML
)
S3method("$", vglm)
@@ -29,6 +44,8 @@ S3method("$", summary.vglm)
S3method("$<-", vglm)
S3method("$<-", summary.vglm)
S3method(formula, vglm)
+S3method(model.matrix, multiple)
+S3method(model.frame, multiple)
S3method(names, relogit)
S3method(names, summary.zelig.relogit)
S3method(names, vglm)
@@ -67,6 +84,7 @@ S3method(print, summary.glm.robust)
S3method(print, summary.MCMCZelig)
S3method(print, summary.zelig)
S3method(print, zelig)
+S3method(repl, default)
S3method(repl, zelig)
S3method(sim, default)
S3method(sim, cond)
@@ -93,10 +111,11 @@ S3method(summary, zelig)
S3method(summary, lm.robust)
S3method(summary, glm.robust)
S3method(terms, vglm)
+S3method(terms, multiple)
S3method(vcov, BetaReg)
S3method(vcov, lm)
S3method(vcov, relogit)
S3method(vcov, lm.robust)
S3method(vcov, glm.robust)
-
+S3method(zelig, default)
diff --git a/R/MCMC/class.ind.R b/R/MCMC/class.ind.R
deleted file mode 100644
index 9c11983..0000000
--- a/R/MCMC/class.ind.R
+++ /dev/null
@@ -1,11 +0,0 @@
-#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/param.MCMCZelig.R b/R/MCMC/param.MCMCZelig.R
deleted file mode 100644
index 007ee15..0000000
--- a/R/MCMC/param.MCMCZelig.R
+++ /dev/null
@@ -1,14 +0,0 @@
-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/MCMC/print.summary.MCMCZelig.R b/R/MCMC/print.summary.MCMCZelig.R
deleted file mode 100644
index d9f634d..0000000
--- a/R/MCMC/print.summary.MCMCZelig.R
+++ /dev/null
@@ -1,12 +0,0 @@
-print.summary.MCMCZelig <- function(object, digits=max(3, getOption("digits") - 3), ...) {
- cat("\nCall: ")
- print(object$call)
- cat("\n", "Iterations = ", object$start, ":", object$end, "\n", sep = "")
- cat("Thinning interval =", object$thin, "\n")
- cat("Number of chains =", object$nchain, "\n")
- cat("Sample size per chain =", (object$end -
- object$start)/object$thin + 1, "\n")
- cat("\n", "Mean, standard deviation, and quantiles for marginal posterior distributions.", "\n")
- print.matrix(round(object$summary, digits=digits))
- cat("\n")
-}
diff --git a/R/MCMC/qi.MCMCZelig.R b/R/MCMC/qi.MCMCZelig.R
deleted file mode 100644
index 924d43c..0000000
--- a/R/MCMC/qi.MCMCZelig.R
+++ /dev/null
@@ -1,393 +0,0 @@
-
-qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
-
- model <- object$zelig
- qi <- list()
- check <- FALSE
-
- if (model %in% c("MCMClogit", "MCMCprobit", "MCMCoprobit", "MCMCmnl"))
- check <- TRUE
-
- if (model %in% c("MCMClogit","MCMCprobit", "MCMCregress",
- "MCMCpoisson","MCMCtobit")) {
-
- if (model == "MCMClogit") {
- 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 == "MCMCprobit") {
- 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 =="MCMCregress") {
- 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 =="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") {
- 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 == "MCMClogit") {
- 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 == "MCMCprobit") {
- 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 == "MCMCregress") {
- 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 == "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
-
- 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("MCMClogit", "MCMCprobit", "MCMCpoisson"))
- {
- 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 =="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
- {
- 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)[[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("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/summary.MCMCZelig.R b/R/MCMC/summary.MCMCZelig.R
deleted file mode 100644
index c297013..0000000
--- a/R/MCMC/summary.MCMCZelig.R
+++ /dev/null
@@ -1,19 +0,0 @@
-
-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/MCMC/zelig2MCMC.R b/R/MCMC/zelig2MCMC.R
deleted file mode 100644
index 134cab3..0000000
--- a/R/MCMC/zelig2MCMC.R
+++ /dev/null
@@ -1,348 +0,0 @@
-zelig2MCMChierEI <- 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)
-}
-
-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)
- 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)
-}
-
-zelig2MCMCprobit <- 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)
-}
-
-zelig2MCMCregress <- 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)
-}
-
-zelig2MCMCpoisson <- 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)
-}
-
-zelig2MCMCtobit <- 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)
-}
-
-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
- {
- 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)
-}
-
-zelig2MCMCoprobit <- 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)
-}
-
-
-
-
-
-
-
-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
- {
- 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)
-}
-
-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
- {
- 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)
-}
-
-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
- {
- 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
deleted file mode 100644
index b345df5..0000000
--- a/R/MCMC/zelig3MCMC.R
+++ /dev/null
@@ -1,135 +0,0 @@
- zelig3MCMCdynamicEI <- zelig3MCMChierEI <- 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
-}
-
-zelig3MCMClogit <- zelig3MCMCoprobit <- zelig3MCMCpoisson <-
- zelig3MCMCmnl <- zelig3MCMCregress <- 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
-}
-
-zelig3MCMCprobit <- 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
-
- }
-
- 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) {
-
- 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
-}
-
-
- 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/MULTIPLE/cmsystemfit.R b/R/MULTIPLE/cmsystemfit.R
deleted file mode 100644
index 09b5de3..0000000
--- a/R/MULTIPLE/cmsystemfit.R
+++ /dev/null
@@ -1,19 +0,0 @@
-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/MULTIPLE/model.matrix.multiple.R b/R/MULTIPLE/model.matrix.multiple.R
deleted file mode 100644
index 5752a5d..0000000
--- a/R/MULTIPLE/model.matrix.multiple.R
+++ /dev/null
@@ -1,67 +0,0 @@
-model.matrix.multiple <- function (object,data,eqn=NULL,...){
-
- # print("model.matrix.multiple is called")
- if(class(object)[[1]]=="formula"){
- terms <-terms(object)
- obj<-"formula"
- }
- else
- {
- obj<-"terms"
- terms<-object
- }
- att<-attributes(terms)
- expVar<-att$term.labels
- nrExpVariables<- length(expVar)
- nrEquations<-length(all.vars((att$variables)[[2]],unique=FALSE))
-
- terms1<-terms
- class(terms1)<-class(terms1)[class(terms1)!="multiple"]
- multiple<-model.matrix.default(terms1,data)
-
-
- attrList<-attributes(multiple)
- if(obj=="formula"){
- if (hasArg(omit))
- omitLst=omit
- else
- omitLst=NULL
- if(hasArg(constrain))
- constrainLst=constrain
- else
- constrainLst=list()
- omitconsLst<-omitconstrain(object,omitLst,constrainLst)
- attrList<-c(attrList,omitconsLst)
- } # end of "if obj="formula""
- else
- {
- if ("omit" %in% names(attributes(object))){
- attrList[["omit"]]<-attributes(object)$omit
- omitAttr<-attributes(object)$omit
- nrEquations<-dim(omitAttr)[[1]]
- }
-
- if ("constrain" %in% names(attributes(object)))
- attrList[["constrain"]]<-attributes(object)$constrain
- }
- if (!is.null(eqn))
- {
-
- eqnattr<-multiple[,c(1:dim(multiple)[[2]])* c(1,as.numeric(!(attrList[["omit"]][eqn,])))]
- attrname<-paste("eqn",eqn,sep="")
- multiple<-eqnattr
- }
- else
- {
- for (i in 1:nrEquations){
- eqntmp<-paste("eqn",i,sep="")
- attrtmp<-c(1:dim(multiple)[[2]])* c(1,as.numeric(!(attrList[["omit"]][i,])))
- attrtmp1<-attrtmp[attrtmp !=0]
- attr(multiple,eqntmp)<-attrtmp1
- }
- }
-
- multiple
-
-}
-
diff --git a/R/MULTIPLE/terms.multiple.R b/R/MULTIPLE/terms.multiple.R
deleted file mode 100644
index 9f719a4..0000000
--- a/R/MULTIPLE/terms.multiple.R
+++ /dev/null
@@ -1,19 +0,0 @@
-terms.multiple <- function (object,omit=NULL, constrain=NULL){
-
- if (any(class(object)=="multiple")){
- terms<-object$terms
- class(terms)<-c(class(terms),"multiple")
- return (terms)
- }
-
- terms<-terms.formula(object)
-
- attrList<-attributes(terms)
- tmp<-omitconstrain(object,omit,constrain)
- attrList[["omit"]]<-tmp[["omit"]]
- attrList[["constrain"]]<-tmp[["constrain"]]
- attributes(terms)<-attrList
- class(terms)<-c(class(terms),"multiple")
- terms
-}
-
diff --git a/R/MULTIPLE/zelig2w2sls.R b/R/MULTIPLE/zelig2w2sls.R
deleted file mode 100644
index 438f507..0000000
--- a/R/MULTIPLE/zelig2w2sls.R
+++ /dev/null
@@ -1,15 +0,0 @@
-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/MULTIPLE_NEW/model.frame.list.R b/R/MULTIPLE_NEW/model.frame.list.R
deleted file mode 100644
index 2782c94..0000000
--- a/R/MULTIPLE_NEW/model.frame.list.R
+++ /dev/null
@@ -1,30 +0,0 @@
-model.frame.list <- function (object,data,...){
- if(class(object)[[1]]=="terms"){
- terms <-object
- }else{
- terms<-terms(object)
- }
- Xnames<-attr(terms,"indVars")
- rhs<-toBuildFormula(Xnames)
- if(!(is.null(rhs)))
- rhs<-(paste("~",rhs))
- else
- stop("NULL dataframe")
- Ynames<-attr(terms,"depVars")
- if(length(Ynames)>1){
- lhs<-toBuildFormula(Ynames,",")
- if (!(is.null(lhs))){
- lhs<-paste("cbind(",lhs)
- lhs<-paste(lhs,")")
- }
- }else{
- lhs=Ynames
- }
- lhs<-as.formula(paste(lhs,rhs))
- Y<-model.frame.default(lhs,data=data)
- result=Y
- # attr(terms,"response")<-1
- attr(result,"terms")<-terms
- result
-}
-
diff --git a/R/MULTIPLE_NEW/model.matrix.list.R b/R/MULTIPLE_NEW/model.matrix.list.R
deleted file mode 100644
index ede51c9..0000000
--- a/R/MULTIPLE_NEW/model.matrix.list.R
+++ /dev/null
@@ -1,171 +0,0 @@
-model.matrix.list <- function (object,data,shape="compact",eqn=NULL,...){
-
- if((shape != "compact") && (shape != "array") && (shape !="stacked"))
- stop("wrong shape argument!")
-
- if(class(object)[[1]]=="terms"){
- terms <-object
- }
- else
- {
- terms<-terms(object)
- }
- if (!(all(eqn %in% names(terms))))
- stop("wrong eqn argument")
-
- constraints<-attr(terms,"constraints")
- intercAttr<-attr(terms,"intercept")
- termlabels<-attr(terms,"term.labels")
- if (is.null(eqn))
- eqn=names(terms)
-
- if (length(eqn)==1)
- shape="compact"
- Xnames<-unique(unlist(termlabels[eqn]))
- rhs<-toBuildFormula(Xnames)
- if(!(is.null(rhs)))
- rhs<-as.formula(paste("~",rhs))
-
- if (shape=="compact"){
- result=model.matrix.default(rhs,data=data)
- if(all(intercAttr==0)){
- result= result[,colnames(result)!="(Intercept)"]
- }
- attr(terms,"parameters")<-colnames(result)
- attr(result,"terms")<-terms
- return(result)
- }
- intercAttr<-intercAttr[eqn]
- termlabels<-termlabels[eqn]
- nrEquations<-length(eqn)
- result<-list()
- for(i in 1:nrEquations)
- result[[eqn[[i]]]]<-c()
-
- ronames<-rownames(data)
- ronr<-nrow(data)
-
- parameters<-c()
- # start with intercept matrix
- matrixInterc<-list()
- howmanyinterc<-length(intercAttr[intercAttr!=0])
-
- if(howmanyinterc>0){
- for (i in 1:nrEquations){
- if(length(termlabels[[eqn[[i]]]])>0){
- matrixInterc[[eqn[[i]]]]<-matrix(0, nrow=ronr,ncol=howmanyinterc,dimnames=list(ronames,names(intercAttr[intercAttr!=0])))
- if(eqn[[i]] %in% names(intercAttr[intercAttr!=0]))
- matrixInterc[[eqn[[i]]]][,eqn[[i]]]<-1
- colnames(matrixInterc[[i]])<- paste("(Intercept)",colnames(matrixInterc[[i]]),sep=":")
- }
- }
- result<-matrixInterc
- parameters<-c(parameters,colnames(matrixInterc[[1]]))
- }
- matrixTmp<-list()
- newXnames<-list()
-
- for(i in 1:nrEquations){
- fi<-toBuildFormula(termlabels[[eqn[[i]]]])
- if(!(is.null(fi))){
- fi<-as.formula(paste("~",fi))
- tmp<-model.matrix.default(fi,data=data)
- matrixTmp[[eqn[[i]]]]<-tmp[,2:ncol(tmp)]
- newXnames[[eqn[[i]]]]<-colnames(matrixTmp[[eqn[[i]]]])
- }
- }
-
- parol<-c()
- if (is.matrix(constraints)){
- constraints<-matrix(constraints[eqn,],nrow=length(eqn),ncol=ncol(constraints),dimnames=list(eqn,colnames(constraints)))
- matrixConstr<-list()
- paramConstr<-c()
-
- for (i in 1:ncol(constraints)){
- if (any(!(is.na(constraints[,i])))){
- paramConstr<-c(paramConstr,colnames(constraints)[[i]])
- }
- }
- for (k in 1:length(paramConstr)){
- tmp<-paramConstr[[k]]
- for (m in 1:nrEquations){
- if(!(is.na(constraints[eqn[[m]],paramConstr[[k]]])))
- tmp<-paste(tmp,eqn[[m]],sep=":")
- }
- parol<-c(parol,tmp)
- }
-
- for (i in 1:nrEquations){
- if(length(termlabels[[eqn[[i]]]])>0){
- matrixConstr[[eqn[[i]]]]<-matrix(0,nrow=ronr,ncol=length(paramConstr),dimnames=list(ronames,paramConstr))
-
- for(j in 1:length(paramConstr)){
- tmpconstrv<-constraints[eqn[[i]],paramConstr[[j]]]
- if(!(is.na(tmpconstrv))){
- matrixConstr[[eqn[[i]]]][,paramConstr[[j]]]<-matrixTmp[[eqn[[i]]]][,tmpconstrv]
- tmp<- newXnames[[eqn[[i]]]]
- newXnames[[eqn[[i]]]]<-tmp[tmp!=tmpconstrv]
- tmp<-matrixTmp[[eqn[[i]]]]
- matrixTmp[[eqn[[i]]]]<-as.matrix(tmp[,newXnames[[eqn[[i]]]]])
- }
- }
- }
- }
- # if there is been not intercept then this would be the first .. otherwise cbind it
- if(length(result)==0)
- result<-matrixConstr
- else
- for(i in 1:nrEquations)
- result[[eqn[[i]]]]<-cbind(result[[eqn[[i]]]],matrixConstr[[eqn[[i]]]])
- parameters<-c(parameters,parol)
- }
-
- for(i in 1:nrEquations){
- tmp<- newXnames[[eqn[[i]]]]
- if(length(tmp)>0){
- newXnames[[eqn[[i]]]]<-paste(tmp,eqn[[i]],sep=":")
- colnames(matrixTmp[[eqn[[i]]]])<-newXnames[[eqn[[i]]]]
- }
- }
- matrixTheRest<-list()
- namesTheRest<-unique(unlist(newXnames ))
-
- for (i in 1:nrEquations){
- if(length(newXnames[[eqn[[i]]]])>0){
- matrixTheRest[[eqn[[i]]]]<-matrix(0,nrow=ronr,ncol=length(namesTheRest),dimnames=list(ronames,namesTheRest))
- matrixTheRest[[eqn[[i]]]][,newXnames[[eqn[[i]]]]]<-matrixTmp[[eqn[[i]]]]
- }
- }
- # if there is been not intercept then this would be the first .. otherwise cbind it
- if(length(result)==0)
- result<-matrixTheRest
- else
- for(i in 1:nrEquations)
- result[[eqn[[i]]]]<-cbind(result[[eqn[[i]]]],matrixTheRest[[eqn[[i]]]])
- parameters<-c(parameters,colnames(matrixTheRest[[1]]))
-
- if(shape=="array"){
- res<-array(0,dim=c(ronr,length(parameters),length(result)),dinames<-list(ronames,colnames(result[[1]]),names(result)))
- for (i in 1:length(result))
- res[,,names(result)[[i]]]<-result[[i]]
- }
-
- if(shape=="stacked"){
- res<-result[[1]]
- if(length(result)>1)
- for(i in 2:length(result))
- res<-rbind(res,result[[i]])
- }
-
- for(i in 1:nrEquations)
- if(length(termlabels[[i]])==0)
- parameters<-c(parameters,eqn[[i]])
-
- attr(terms,"parameters")<-parameters
-
-
-
- attr(res,"terms")<-terms
- return(res)
-}
-
diff --git a/R/MULTIPLE_NEW/parse.formula.R b/R/MULTIPLE_NEW/parse.formula.R
deleted file mode 100644
index 1f21f48..0000000
--- a/R/MULTIPLE_NEW/parse.formula.R
+++ /dev/null
@@ -1,90 +0,0 @@
-parse.formula<-function( listofformulas, req=NULL, opt=NULL, fixed=NULL){
-
-nrrho<-0 # number of opt ???? replace nrrho with nopts
-nreq<-0 # final number of equations (for req) i.e. rep contain more then one eq)
-eqns<-list() # list of req equations
-rho<-list() # list of opt equations ??? replace rho with opt???
-fix<-list()
-res<-list() # the final result : list of equation + attributes
-
-#if(!(is.null(fixed)) && fixed=="rho")
-# stop("Please chose a different value for fixed parameter")
-
-if(class(listofformulas)=="formula")
- listofformulas<-list(listofformulas)
-
-
-nreqns <-length(listofformulas) # total number of equations
-for (i in 1:nreqns){
- eqni<-listofformulas[[i]]
- if (length(eqni)==3){ # on left hand side we have something
- lhs<-deparse(eqni[[2]])
- rhs<-deparse(eqni[[3]])
- if(substr(lhs,1,5)!="cbind"){
- nreq=nreq+1
- e<-paste(lhs,"~",sep="")
- eqns[[nreq]]<-as.formula(paste(e,rhs,sep=""))
- }
- else
- {
- lhs<-eqni[[2]]
- g<- all.vars(lhs) # how many eq we have inside cbind
- for (j in 1:length(g)){
- nreq=nreq+1
- e<-paste(g[[j]],"~",sep="")
- eqns[[nreq]]<-as.formula(paste(e,rhs,sep=""))
- }
- }
- }
- else # if we dont have anyting means is rho
- {
- rhs<-deparse(eqni[[2]])
- nrrho=nrrho+1
- rho[[nrrho]]<-as.formula(paste("~",rhs,sep=""))
- }
-}
-
-if (is.null(opt)){
- if (nrrho !=0)
- stop("number of equations does not match model inputs!")
-}
-else
- {
- if (nrrho < length(opt))
- if (length(opt)==1){
- nrrho=nrrho+1
- rho[[nrrho]]<-as.formula("~1")
- }
- else
- stop("number of equations does not match model inputs!")
- names(rho)<-opt
- if (nrrho > length(opt))
- if(length(opt)==1)
- names(rho)<-paste(opt,1:nrrho,sep="")
- else
- stop("number of equations does not match model inputs!")
- }
-
-if(!(is.null(fixed)))
-for(i in 1:length(fixed)){
- fix[[fixed[[i]]]]<-as.formula("~1")
-}
-
-if (is.null(req)){
- if (nreq !=0)
- stop("number of equations does not match model inputs!")
-}
-else
- {
- if (nreq < length(req))
- stop("number of equations does not match model inputs!")
- names(eqns)<-req
- if (nreq > length(req))
- if(length(req)==1)
- names(eqns)<-paste(req,1:nreq,sep="")
- else
- stop("number of equations does not match model inputs!")
- }
-res<-c(eqns,rho,fix)
-return(res)
-}
diff --git a/R/MULTIPLE_NEW/parse.par.R b/R/MULTIPLE_NEW/parse.par.R
deleted file mode 100644
index 3bb1a9c..0000000
--- a/R/MULTIPLE_NEW/parse.par.R
+++ /dev/null
@@ -1,74 +0,0 @@
-parse.par <- function(par, terms, eqn, shape = NULL) {
- "%w/o%" <- function(x,y) x[!x %in% y]
- if (is.null(shape)) {
- if (any(class(terms) == "list"))
- shape <- "matrix"
- else
- shape <- "vector"
- }
- if (!shape %in% c("matrix", "vector"))
- stop("not a valid 'shape' for parameters. Choose from \"matrix\" or \"vector\".")
- if (any(class(terms) == "list")) {
- idx <- make.parameters(terms = terms, shape = "vector")
- mat <- t(make.parameters(terms = terms, shape = "matrix"))
- syst <- rownames(mat)
- if (length(syst) == 1)
- shape <- "vector"
- par.names <- unique(na.omit(c(mat)))
- ancil <- idx %w/o% par.names
- if (any(eqn %in% ancil)) {
- if (any(eqn %in% syst)) {
- stop(" eqn cannot include both systematic and ancillary \n parameters at the same time.")
- }
- else
- ret.val <- par[idx %in% eqn]
- }
- else { ## if eqn to be returned is a systematic component
- if (length(eqn) > 1) {
- sys.idx <- sort(pmatch(eqn, syst))
- var.idx <- sort(pmatch(eqn, idx))
- }
- else {
- sys.idx <- sort(grep(eqn, syst))
- var.idx <- sort(grep(eqn, idx))
- }
- subs <- mat[sys.idx, , drop = FALSE]
- out <- matrix(0, nrow = nrow(subs), ncol = ncol(subs),
- dimnames = dimnames(subs))
- nas <- which(is.na(subs), arr.ind = TRUE)
- not.na <- which(!is.na(subs), arr.ind = TRUE)
- const <- attr(terms, "constraints")
- if (is.logical(const)) { ## for no constraints
- out[not.na] <- par[var.idx]
- if (shape == "matrix")
- ret.val <- t(out)
- else {
- out[nas] <- NA
- ret.val <- na.omit(c(t(out)))
- }
- }
- else { ## if constraints
- for (i in 1:ncol(const)) {
- idx <- idx %w/o% ancil
- const.loc <- which(subs == colnames(const)[i], arr.ind = TRUE)
- for (j in 1:nrow(const.loc)) {
- rm.idx <- apply(not.na, 1, identical, const.loc[j,])
- not.na[rm.idx,] <- NA
- }
- out[const.loc] <- par[idx == colnames(const)[i]]
- }
- not.na <- na.omit(not.na)
- out[not.na] <- par[1:nrow(not.na)]
- if (shape == "matrix")
- ret.val <- t(out)
- else {
- out[nas] <- NA
- ret.val <- unique(na.omit(c(t(out)))) ## FIX
- }
- }
- }
- }
- ret.val
-}
-
-
diff --git a/R/MULTIPLE_NEW/tag.R b/R/MULTIPLE_NEW/tag.R
deleted file mode 100644
index 1f88071..0000000
--- a/R/MULTIPLE_NEW/tag.R
+++ /dev/null
@@ -1,5 +0,0 @@
-tag <- function (x,y){
-return(x)
-
-}
-
diff --git a/R/MULTIPLE_NEW/terms.list.R b/R/MULTIPLE_NEW/terms.list.R
deleted file mode 100644
index 46abcb7..0000000
--- a/R/MULTIPLE_NEW/terms.list.R
+++ /dev/null
@@ -1,113 +0,0 @@
-terms.list<-function(object,...){
-
-nreq<- 0
-reqEqns<-list()
-optEqns<-list()
-constr<-list()
-XconstrEqn<-list()
-variables<-list()
-termlabels<-list()
-depVars<-c()
-nrConstr<-0
-namesConstr<-c()
-namesOfEquations<- names(object)
-nrEquations <-length(object)
-nrEquationsNew<-0
-objectNew<-list()
-intercAttr<-list()
-
-trim<-function(s){
- return (sub(" *([^ ]+) *", "\\1", s))
-}
-
-
-parseTag<-function(s){
- rml<-sub(".*\\(","",s)
- rmr<-sub("\\).*","",rml)
- r<-unlist(strsplit(rmr,",",fixed=TRUE))
- r<-trim(r)
- if(substr(r[[2]],1,1)!="\"" || substr(r[[2]],nchar(r[[2]]),nchar(r[[2]]))!="\"")
- stop( "the name of the parameters should be label (enclosed in quotation marks ...)")
- r[[2]]<- substr(r[[2]],2,nchar(r[[2]])-1)
- return(r)
-}
-
-
-
-for (i in 1:nrEquations){
- TT<-terms.formula(object[[i]], specials="tag")
- eqni<-object[[i]]
- namei<-namesOfEquations[[i]]
- tagattr<-attr(TT,"specials")$tag
- hastag<-!(is.null(tagattr))
- if (hastag){
- constrTmp<-c()
- for(j in 1:length(tagattr)){
- nrConstr<-nrConstr+1
- if(length(eqni)==3)
- ind<-tagattr[[j]]-1
- else
- ind<-tagattr[[j]]
- parsedTag<-parseTag(attr(TT,"term.labels")[[ind]])
- namesConstr<-c(namesConstr,parsedTag[[2]])
- constr[[nrConstr]]<-c(i,parsedTag[[2]],parsedTag[[1]])
- attr(TT,"term.labels")[[ind]]<-parsedTag[[1]]
- ind<-tagattr[[j]]+1
- attr(TT,"variables")[[ind]]<-as.name(parsedTag[[1]])
- constrTmp<-c(constrTmp,parsedTag[[1]])
- }
- XconstrEqn[[i]]<-constrTmp
- }
-
- if (length(eqni)==3){
- nrEquationsNew<-nrEquationsNew+1
- objectNew[[namei]]<-eqni
- nreq=nreq+1
- reqEqns[[namei]]<-eqni
- depVars<-c(depVars,deparse(eqni[[2]]))
- }else{
- nrEquationsNew<-nrEquationsNew+1
- objectNew[[namei]]<-eqni
- optEqns[[namei]]<-eqni
- }
- variables[[namei]]<-attr(TT,"variables")
- termlabels[[namei]]<-attr(TT,"term.labels")
- intercAttr[[namei]]<-attr(TT,"intercept")
-}
-
-namesOfEquations<-names(objectNew)
-
-myattr<-list()
-result<-objectNew
-reqnames<-names(reqEqns)
-
-if(length(constr)>0){
- namesConstr<-unique(namesConstr)
- constraints<-matrix(NA,nrow=nrEquationsNew,ncol=length(namesConstr),dimnames=list(namesOfEquations,namesConstr))
- for(i in 1:length(constr)){
- constri<-constr[[i]]
- eqind<-constri[[1]]
- eq<-namesOfEquations[as.numeric(eqind)]
- lab<-constri[[2]]
- constraints[eq,lab]<-constri[[3]]
- }
-}else
-constraints<-FALSE
-
-indVars<-unique(unlist(termlabels))
-
-myattr$reqEqns<-reqEqns
-myattr$optEqns<-optEqns
-myattr$variables<-variables
-myattr$term.labels<-termlabels
-myattr$indVars<-indVars
-
-myattr$depVars<-unlist(depVars)
-myattr$constraints<-constraints
-myattr$response<-1
-myattr$intercept<-intercAttr
-attributes(result)<-myattr
-names(result)<-namesOfEquations
-class(result)<-c("terms","list")
-return(result)
-}
diff --git a/R/MULTIPLE_NEW/toBuildFormula.R b/R/MULTIPLE_NEW/toBuildFormula.R
deleted file mode 100644
index fb815d4..0000000
--- a/R/MULTIPLE_NEW/toBuildFormula.R
+++ /dev/null
@@ -1,17 +0,0 @@
-toBuildFormula<-function(Xnames,sepp="+"){
- lng<-length(Xnames)
- rhs=NULL
- if (lng!=0){
- if(lng==1){
- rhs=Xnames
- }else{
- for (j in 1:(lng-1)){
- rhs<-paste(rhs,as.name(Xnames[[j]]))
- rhs<-paste(rhs,sepp)
- }
- rhs<-paste(rhs,Xnames[[lng]])
- }
- }
-
- return (rhs)
-}
diff --git a/R/categories.R b/R/categories.R
new file mode 100644
index 0000000..b1f80f0
--- /dev/null
+++ b/R/categories.R
@@ -0,0 +1,9 @@
+categories <-function(){
+list(continuous="Regression Models for Continous Dependent Variables",
+ dichotomous="Models for Dichotomous Dependent Variables",
+ ordinal="Models for Ordinal Dependent Variables",
+ bivariate.dichotomous="Models for pairs of Dichotomous Dependent Variables",
+ multinomial="Multinomial Choice Models",
+ event.count="Event Count Models",
+ censored="Models for Censored (and Duration) Variables")
+}
diff --git a/R/check.describe.R b/R/check.describe.R
new file mode 100644
index 0000000..c612263
--- /dev/null
+++ b/R/check.describe.R
@@ -0,0 +1,50 @@
+check.describe<-function(mymodel){
+ firstLvlNames<-c("category","description","package","parameters")
+ fn <- paste("describe.", mymodel, sep = "")
+ if (!exists(fn))
+ stop("Nothing to check ... The function describe.",mymodel,"does not exist")
+ z<-do.call(fn,list())
+
+ #any extra name in the list??
+ whiche<-which(!(names(z) %in% firstLvlNames))
+ if (length(whiche)!=0){
+ tmpstr<-names(z)[[whiche[[1]]]]
+ if(length(whiche)>1)
+ for(i in 2:length(whiche))
+ tmpstr<-paste(tmpstr,names(z)[[whiche[[i]]]],sep=",")
+ stop ("Unknown names in your list: ",tmpstr)
+ }
+ errmsg<-" is missing. It's required ..."
+ if(is.null(z$category))
+ stop("\"category\"", errmsg)
+ else
+ if(!(z$category %in% names(categories())))
+ stop("unknown category \"",z$category, "\"")
+ if(is.null(z$parameters)) stop("\"parameters\"",errmsg)
+
+ for (i in length(z$parameters)){
+ eqns<-z$parameters[[i]]$equations
+ if(is.null(eqns)) stop("\"equations\"",errmsg)
+ if(length(eqns)!=2) stop("equations must be an vector of length 2")
+ if(!(eqns[[2]] <= 999 || !(is.finite(eqns[[2]]))) ) stop("The maximum number of equations for each paramter should be <=999 or \"Inf\"..")
+
+ tags<-z$parameters[[i]]$tagsAllowed
+ if (is.null(tags)) stop ("\"tagsAllowed\"",errmsg)
+ if(!is.logical(tags)) stop("\"tagsAllowed\" must have a logical value (\"TRUE\" or \"FALSE\")")
+
+ dep<-z$parameters[[i]]$depVar
+ if (is.null(dep)) stop ("\"depVar\"",errmsg)
+ if(!is.logical(dep)) stop("\"depVar\" must have a logical value (\"TRUE\" or \"FALSE\")")
+
+ exp<-z$parameters[[i]]$expVar
+ if (is.null(exp)) stop ("\"expVar\"",errmsg)
+ if(!is.logical(exp)) stop("\"expVar\" must have a logical value (\"TRUE\" or \"FALSE\")")
+
+ tags<-z$parameters[[i]]$tagsAllowed
+ if (is.null(tags)) stop ("\"tagsAllowed\"",errmsg)
+ if(!is.logical(tags)) stop("\"tagsAllowed\" must have a logical value (\"TRUE\" or \"FALSE\")")
+ }
+
+ cat("Congratulation, your function \"",fn, "\" passed this test\n")
+}
+
diff --git a/R/cmvglm.R b/R/cmvglm.R
index dfb564c..807594f 100644
--- a/R/cmvglm.R
+++ b/R/cmvglm.R
@@ -1,20 +1,36 @@
-cmvglm <- function(formula, model, constrain, omit, constant, ndim){
+cmvglm <- function(formula, model, constant, ndim,data=NULL, fact=NULL){
+
+ toBuildFormula<-function(Xnames,sepp="+"){
+ lng<-length(Xnames)
+ rhs<-NULL
+ if (lng!=0){
+ if(lng==1){
+ rhs=Xnames
+ }else{
+ for (j in 1:(lng-1)){
+ rhs<-paste(rhs,as.name(Xnames[[j]]))
+ rhs<-paste(rhs,sepp)
+ }
+ rhs<-paste(rhs,Xnames[[lng]])
+ }
+ }
+ return (rhs)
+ }
tt<-terms(formula)
- vars<-c("(Intercept)", attr(tt, "term.labels"))
+ attr(tt,"systEqns")<-names(formula)
+ p<-make.parameters(tt,shape="matrix")
+ vars<-rownames(p)
cm<-vector("list", length(vars))
names(cm)<-vars
+
for(i in 1:length(cm))
cm[[i]]<-diag(1, ndim)
- if(!is.null(constrain)){
- tmp <- sort(names(constrain))
- for (i in 1:length(constrain[[tmp[1]]])) {
- ci <- NULL
- counter <- c("1", "2", "3")
- for (j in 1:3)
- if (is.null(constrain[[counter[j]]]))
- ci <- c(ci, NA)
- else
- ci <- c(ci, constrain[[counter[j]]][i])
+
+ constrain<-attr(tt,"constraints")
+ if(!is.logical(constrain)){
+ tmp <- sort(colnames(constrain))
+ for (i in 1:length(tmp)) {
+ ci<-constrain[,i]
if (is.null(na.omit(ci)) || length(unique(na.omit(ci)))!=1)
stop("invalid input for constrain")
minj <- match(FALSE, is.na(ci))
@@ -26,14 +42,13 @@ cmvglm <- function(formula, model, constrain, omit, constant, ndim){
}
}
}
- if(!is.null(omit)){
- tmp <- sort(names(omit))
- for(i in 1:length(omit)){
- omiti <- omit[[tmp[i]]]
- for (j in 1:length(omiti))
- cm[[pmatch(omiti[j], names(cm))]][i,i]<-0
+ for(i in rownames(p)){
+ for(j in 1:ncol(p)){
+ if(is.na(p[i,j]))
+ cm[[i]][j,j]<-0
}
}
+
if(!is.null(constant))
for(i in 1:length(constant))
for(j in 1:length(cm))
@@ -41,5 +56,25 @@ cmvglm <- function(formula, model, constrain, omit, constant, ndim){
cm[[j]][constant[i],]<-matrix(0, ncol=ncol(cm[[j]]))
for(i in 1:length(cm))
cm[[i]]<-as.matrix(cm[[i]][,apply(cm[[i]], 2, sum)!=0])
+ rhs<-toBuildFormula(attr(tt,"indVars"))
+ if(!(is.null(rhs)))
+ rhs<-(paste("~",rhs))
+ else
+ rhs<-"~1"
+ Ynames<-unlist(attr(tt,"depVars"))
+ if(!is.null(fact))
+ lhs<-fact
+ else{
+ if(length(Ynames)>1){
+ lhs<-toBuildFormula(Ynames,",")
+ if (!(is.null(lhs))){
+ lhs<-paste("cbind(",lhs)
+ lhs<-paste(lhs,")")
+ }
+ }else{
+ lhs=Ynames
+ }
+ }
+ formula<-as.formula(paste(lhs,rhs))
list("formula"=formula, "constraints"=cm)
}
diff --git a/R/current.packages.R b/R/current.packages.R
new file mode 100644
index 0000000..e820aaa
--- /dev/null
+++ b/R/current.packages.R
@@ -0,0 +1,51 @@
+current.packages <- function(package){
+
+ required.packages <- function(pack) {
+ mylib <- dirname(system.file(package = pack))
+ description <- packageDescription(pack, lib = mylib)
+ depends <- description$Depends
+ if (!is.null(depends)) {
+ depends <- strsplit(depends, ", ")[[1]]
+ Rdepends <- pmatch("R (", depends)
+ if (is.na(Rdepends))
+ Rdepends <- pmatch("R(", depends)
+ if (!is.na(Rdepends))
+ depends <- depends[-pmatch("R (", depends)]
+ }
+ suggests <- description$Suggests
+ if (!is.null(suggests))
+ suggests <- strsplit(suggests, ", ")[[1]]
+ total <- c(depends, suggests)
+ if (!is.null(total)) {
+ conditions <- grep(")", total)
+ if (length(conditions) > 0) {
+ for (i in conditions)
+ total[i] <- strsplit(total[i], " \\(")[[1]][1]
+ }
+ return(total)
+ }
+ else
+ return(NULL)
+ }
+ old <- packages <- required.packages(package)
+
+ for (zpack in packages) {
+ new <- required.packages(zpack)
+ tmp <- new[!(new %in% old)]
+ old <- packages <- c(packages, tmp)
+ }
+
+ ver <- array()
+ for (zpack in na.omit(packages)) {
+ mylib <- dirname(system.file(package = zpack))
+ ver[zpack] <- packageDescription(zpack, lib = mylib)$Ver
+ }
+ ver[1] <- paste(paste(paste(R.Version()$major, R.Version()$minor, sep = "."),
+ R.Version()$status, sep = " "),
+ R.Version()$svn, sep = " svn: ")
+ names(ver)[1] <- "R"
+ vv <- as.matrix(ver)
+ colnames(vv) <- "Version"
+ noquote(vv)
+}
+
diff --git a/R/describe.blogit.R b/R/describe.blogit.R
new file mode 100644
index 0000000..939424a
--- /dev/null
+++ b/R/describe.blogit.R
@@ -0,0 +1,19 @@
+describe.blogit<-function(){
+category <- "dichotomous"
+description <- "Bivariate Logistic Regression for Dichotomous Dependent Variables"
+package <-list( name ="VGAM",
+ version ="0.6",
+ CRAN ="http://www.stat.auckland.ac.nz/~yee"
+ )
+parameters<-list(mu="mu",phi="phi")
+parameters$mu<-list(equations=c(2,2),
+ tagsAllowed=TRUE,
+ depVar=TRUE,
+ expVar=TRUE)
+
+parameters$phi<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=FALSE,
+ expVar=TRUE)
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.bprobit.R b/R/describe.bprobit.R
new file mode 100644
index 0000000..95a003d
--- /dev/null
+++ b/R/describe.bprobit.R
@@ -0,0 +1,19 @@
+describe.bprobit<-function(){
+category <- "dichotomous"
+description <- "Bivariate Probit Regression for Dichotomous Dependent Variables"
+package <-list( name ="VGAM",
+ version ="0.6",
+ CRAN ="http://www.stat.auckland.ac.nz/~yee"
+ )
+parameters<-list(mu="mu", rho="rho")
+parameters$mu<-list(equations=c(2,2),
+ tagsAllowed=TRUE,
+ depVar=TRUE,
+ expVar=TRUE)
+
+parameters$rho<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=FALSE,
+ expVar=TRUE)
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.default.R b/R/describe.default.R
new file mode 100644
index 0000000..e4cbe9b
--- /dev/null
+++ b/R/describe.default.R
@@ -0,0 +1,12 @@
+describe.default<-function(){
+category <- "Dichotomous"
+description <- "A statistical model"
+
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE)
+
+list(category=category,description=description,parameters=parameters)
+}
diff --git a/R/describe.exp.R b/R/describe.exp.R
new file mode 100644
index 0000000..0421b39
--- /dev/null
+++ b/R/describe.exp.R
@@ -0,0 +1,17 @@
+describe.exp<-function(){
+category <- "censored"
+description <- "Exponential Regression for Duration Dependent Variables"
+package <-list( name ="survival",
+ version ="2.0",
+ CRAN =NA
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE,
+ specialFunction="Surv",
+ varInSpecialFunction=c(2,2))
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.factor.bayes.R b/R/describe.factor.bayes.R
new file mode 100644
index 0000000..2d925ae
--- /dev/null
+++ b/R/describe.factor.bayes.R
@@ -0,0 +1,15 @@
+describe.factor.bayes<-function(){
+category <- "censored"
+description <- "Bayesian Factor Analysis"
+package <-list( name ="MCMCpack",
+ version ="0.6"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE,
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.factor.mix.R b/R/describe.factor.mix.R
new file mode 100644
index 0000000..b605f93
--- /dev/null
+++ b/R/describe.factor.mix.R
@@ -0,0 +1,17 @@
+describe.factor.mix<-function(){
+category <- "censored"
+description <- "Mixed Data Factor Analysis"
+package <-list( name ="MCMCpack",
+ version ="0.6"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=FALSE,
+ specialFunction="cbind",
+ varInSpecialFunction=c(1,Inf)
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.factor.ord.R b/R/describe.factor.ord.R
new file mode 100644
index 0000000..a84885c
--- /dev/null
+++ b/R/describe.factor.ord.R
@@ -0,0 +1,17 @@
+describe.factor.ord<-function(){
+category <- "censored"
+description <- "Ordinal Data Factor Analysis"
+package <-list( name ="MCMCpack",
+ version ="0.6"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=FALSE,
+ specialFunction="cbind",
+ varInSpecialFunction=c(1,Inf)
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.gamma.R b/R/describe.gamma.R
new file mode 100644
index 0000000..b7bb59c
--- /dev/null
+++ b/R/describe.gamma.R
@@ -0,0 +1,15 @@
+describe.gamma<-function(){
+category <- "continuous"
+description <- "Gamma Regression for Continuous, Positive Dependent Variables"
+package <-list( name ="stats",
+ version ="0.1"
+ )
+parameters<-list(lambda="lambda")
+parameters$lambda<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.irt1d.R b/R/describe.irt1d.R
new file mode 100644
index 0000000..c18711d
--- /dev/null
+++ b/R/describe.irt1d.R
@@ -0,0 +1,17 @@
+describe.factor.irt1d<-function(){
+category <- "censored"
+description <- "One Dimensional Item Response Model"
+package <-list( name ="MCMCpack",
+ version ="0.6"
+ )
+parameters<-list(pi="pi")
+parameters$pi<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=FALSE,
+ specialFunction="cbind",
+ varInSpecialFunction=c(1,Inf)
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.irtkd.R b/R/describe.irtkd.R
new file mode 100644
index 0000000..cabfefd
--- /dev/null
+++ b/R/describe.irtkd.R
@@ -0,0 +1,17 @@
+describe.irtkd<-function(){
+category <- "censored"
+description <- "K-Dimensional Item Response Model"
+package <-list( name ="MCMCpack",
+ version ="0.6"
+ )
+parameters<-list(pi="pi")
+parameters$pi<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=FALSE,
+ specialFunction="cbind",
+ varInSpecialFunction=c(2,Inf)
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.logit.R b/R/describe.logit.R
new file mode 100644
index 0000000..9a88102
--- /dev/null
+++ b/R/describe.logit.R
@@ -0,0 +1,15 @@
+describe.logit<-function(){
+category <- "dichotomous"
+description <- "Logistic Regression for Dichotomous Dependent Variables"
+package <-list( name ="stats",
+ version ="0.1"
+ )
+parameters<-list(pi="pi")
+parameters$pi<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.logit.bayes.R b/R/describe.logit.bayes.R
new file mode 100644
index 0000000..a90b3b1
--- /dev/null
+++ b/R/describe.logit.bayes.R
@@ -0,0 +1,15 @@
+describe.logit.bayes<-function(){
+category <- "dichotomous"
+description <- "Bayesian Logistic Regression for Dichotomous Dependent Variables"
+package <-list( name ="MCMCpack",
+ version ="0.6",
+ )
+parameters<-list(pi="pi")
+parameters$pi<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.lognorm.R b/R/describe.lognorm.R
new file mode 100644
index 0000000..57d6e31
--- /dev/null
+++ b/R/describe.lognorm.R
@@ -0,0 +1,17 @@
+describe.lognorm<-function(){
+category <- "censored"
+description <- "Log-Normal Regression for Duration Dependent Variables"
+package <-list( name ="survival",
+ version ="2.2"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE,
+ specialFunction="Surv",
+ varInSpecialFunction=c(2,2)
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.ls.R b/R/describe.ls.R
new file mode 100644
index 0000000..0be84f5
--- /dev/null
+++ b/R/describe.ls.R
@@ -0,0 +1,15 @@
+describe.ls<-function(){
+category <- "continuous"
+description <- "Least Squares Regression for Continuous Dependent Variables"
+package <-list( name ="stats",
+ version ="0.1"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.mlogit.R b/R/describe.mlogit.R
new file mode 100644
index 0000000..1faafcf
--- /dev/null
+++ b/R/describe.mlogit.R
@@ -0,0 +1,18 @@
+describe.mlogit<-function(){
+category <- "multinomial"
+description <- "Multinomial Logistic Regression for Dependent Variables with Unordered Categorical Values"
+package <-list( name ="VGAM",
+ version ="0.6",
+ CRAN="http://www.stat.auckland.ac.nz/~yee"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,Inf),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE,
+ specialFunction="as.factor",
+ varInSpecialFunction=c(1,1)
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.mlogit.bayes.R b/R/describe.mlogit.bayes.R
new file mode 100644
index 0000000..383681c
--- /dev/null
+++ b/R/describe.mlogit.bayes.R
@@ -0,0 +1,17 @@
+describe.mlogit.bayes<-function(){
+category <- "multinomial"
+description <- "Bayesian Multinomial Logistic Regression for Dependent Variables with Unordered Categorical Values"
+package <-list( name ="MCMCpack",
+ version ="0.6"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE,
+ specialFunction="as.factor",
+ varInSpecialFunction=c(1,1)
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.negbin.R b/R/describe.negbin.R
new file mode 100644
index 0000000..b6c72b8
--- /dev/null
+++ b/R/describe.negbin.R
@@ -0,0 +1,15 @@
+describe.negbin<-function(){
+category <- "event.count"
+description <- "Negative Binomial Regression for Event Count Dependent Variables"
+package <-list( name ="MASS",
+ version ="0.1"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.normal.R b/R/describe.normal.R
new file mode 100644
index 0000000..112fcbc
--- /dev/null
+++ b/R/describe.normal.R
@@ -0,0 +1,15 @@
+describe.normal<-function(){
+category <- "countinuous"
+description <- "Normal Regression for Continuous Dependent Variables"
+package <-list( name ="stats",
+ version ="0.1"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.normal.bayes.R b/R/describe.normal.bayes.R
new file mode 100644
index 0000000..e62cdfb
--- /dev/null
+++ b/R/describe.normal.bayes.R
@@ -0,0 +1,15 @@
+describe.normal.bayes<-function(){
+category <- "countinuous"
+description <- "Bayesian Normal Linear Regression"
+package <-list( name ="MCMCpack",
+ version ="0.6"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.ologit.R b/R/describe.ologit.R
new file mode 100644
index 0000000..8e358d9
--- /dev/null
+++ b/R/describe.ologit.R
@@ -0,0 +1,17 @@
+describe.ologit<-function(){
+category <- "ordinal"
+description <- "Ordinal Logistic Regression for Ordered Categorical Dependent Variables"
+package <-list( name ="MASS",
+ version ="0.1"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE,
+ specialFunction="as.factor",
+ varInSpecialFunction=c(1,1)
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.oprobit.R b/R/describe.oprobit.R
new file mode 100644
index 0000000..575dec8
--- /dev/null
+++ b/R/describe.oprobit.R
@@ -0,0 +1,19 @@
+describe.oprobit<-function(){
+
+category <- "ordinal"
+description <- "Ordinal Probit Regression for Ordered Categorical Dependent Variables"
+package <-list( name ="MASS",
+ version ="0.1"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE,
+ specialFunction="as.factor",
+ varInSpecialFunction=c(1,1)
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+
+}
diff --git a/R/describe.oprobit.bayes.R b/R/describe.oprobit.bayes.R
new file mode 100644
index 0000000..e54ba0e
--- /dev/null
+++ b/R/describe.oprobit.bayes.R
@@ -0,0 +1,17 @@
+describe.oprobit.bayes<-function(){
+category <- "ordinal"
+description <- "Bayesian Ordered Probit Regression"
+package <-list( name ="MCMCpack",
+ version ="0.6"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE,
+ specialFunction="as.factor",
+ varInSpecialFunction=c(1,1)
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.poisson.R b/R/describe.poisson.R
new file mode 100644
index 0000000..19e82e3
--- /dev/null
+++ b/R/describe.poisson.R
@@ -0,0 +1,15 @@
+describe.poisson<-function(){
+category <- "event.count"
+description <- "Poisson Regression for Event Count Dependent Variables"
+package <-list( name ="stats",
+ version ="0.1"
+ )
+parameters<-list(lambda="lambda")
+parameters$lambda<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.poisson.bayes.R b/R/describe.poisson.bayes.R
new file mode 100644
index 0000000..d31b69e
--- /dev/null
+++ b/R/describe.poisson.bayes.R
@@ -0,0 +1,15 @@
+describe.poisson.bayes<-function(){
+category <- "event.count"
+description <- "Ordinal Probit Regression for Ordered Categorical Dependent Variables"
+package <-list( name ="MCMCpack",
+ version ="0.6"
+ )
+parameters<-list(lambda="lambda")
+parameters$lambda<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.probit.R b/R/describe.probit.R
new file mode 100644
index 0000000..8dc2631
--- /dev/null
+++ b/R/describe.probit.R
@@ -0,0 +1,15 @@
+describe.probit<-function(){
+category <- "dichotomous"
+description <- "Probit Regression for Dichotomous Dependent Variables"
+package <-list( name ="stats",
+ version ="0.1"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.probit.bayes.R b/R/describe.probit.bayes.R
new file mode 100644
index 0000000..a17dd50
--- /dev/null
+++ b/R/describe.probit.bayes.R
@@ -0,0 +1,15 @@
+describe.probit.bayes<-function(){
+category <- "dichotomous"
+description <- "Bayesian Probit Regression for Dichotomous Dependent Variables"
+package <-list( name ="MCMCpack",
+ version ="0.6"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.relogit.R b/R/describe.relogit.R
new file mode 100644
index 0000000..014cd6a
--- /dev/null
+++ b/R/describe.relogit.R
@@ -0,0 +1,15 @@
+describe.relogit<-function(){
+category <- "dichotomous"
+description <- "Rare Events Logistic Regression for Dichotomous Dependent Variables"
+package <-list( name ="stats",
+ version ="0.1"
+ )
+parameters<-list(pi="pi")
+parameters$pi<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.tobit.bayes.R b/R/describe.tobit.bayes.R
new file mode 100644
index 0000000..469d98f
--- /dev/null
+++ b/R/describe.tobit.bayes.R
@@ -0,0 +1,15 @@
+describe.tobit.bayes<-function(){
+category <- "censored"
+description <- "Bayesian Linear Regression for a Censored Dependent Variable"
+package <-list( name ="MCMCpack",
+ version ="0.6"
+ )
+parameters<-list(mu="mu")
+parameters$mu<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/describe.weibull.R b/R/describe.weibull.R
new file mode 100644
index 0000000..f7a12f6
--- /dev/null
+++ b/R/describe.weibull.R
@@ -0,0 +1,17 @@
+describe.weibull<-function(){
+category <- "censored"
+description <- "Weibull Regression for Duration Dependent Variables"
+package <-list( name ="survival",
+ version ="2.2"
+ )
+parameters<-list(lambda="lambda")
+parameters$lambda<-list(equations=c(1,1),
+ tagsAllowed=FALSE,
+ depVar=TRUE,
+ expVar=TRUE,
+ specialFunction="Surv",
+ varInSpecialFunction=c(2,2)
+ )
+
+list(category=category,description=description,package=package,parameters=parameters)
+}
diff --git a/R/MULTIPLE_NEW/make.parameters.R b/R/make.parameters.R
similarity index 66%
rename from R/MULTIPLE_NEW/make.parameters.R
rename to R/make.parameters.R
index f44349e..29a9924 100644
--- a/R/MULTIPLE_NEW/make.parameters.R
+++ b/R/make.parameters.R
@@ -1,39 +1,41 @@
-make.parameters <- function(terms, shape = "vector", ancillary = TRUE) {
+make.parameters <- function(terms, shape = "vector", ancillary = TRUE,eqns=NULL) {
if (!shape %in% c("matrix", "vector"))
stop("not a valid 'shape' for parameters. Choose from \"matrix\" or \"vector\".")
- ints <- attr(terms, "intercept")
- eqns <- names(ints)
- labs <- attr(terms, "term.labels")
+ #comment
+ if(is.null(eqns))
+ eqns<-names(terms)
+ ints <- attr(terms, "intercept")[eqns]
+ labs <- attr(terms, "term.labels")[eqns]
const <- attr(terms, "constraints")
- const[const == 0] <- NA
for (i in 1:length(eqns)) {
if (ints[[i]] == 1)
labs[[i]] <- c("(Intercept)", labs[[i]])
}
- "%w/o%" <- function(x,y) x[!x %in% y]
- fixed <- which(labs == "(Intercept)")
- syst <- 1:length(eqns) %w/o% fixed
+ fixed<-eqns[eqns %in% attr(terms,"ancilEqns")]
+ syst<-eqns[eqns %in% attr(terms,"systEqns")]
+# syst<-eqns
vars <- unique(unlist(labs))
pars <- matrix(NA, ncol = length(syst), nrow = length(vars))
- colnames(pars) <- eqns[syst]
+ colnames(pars) <- syst
rownames(pars) <- vars
for (i in syst) {
idx <- which(!is.na(match(vars, labs[[i]])))
- pars[idx,i] <- paste(labs[[i]], eqns[i], sep = ":")
+ pars[idx,i] <- paste(labs[[i]], i, sep = ":")
}
if (!is.logical(const)) {
+ const <- attr(terms, "constraints")[syst,,drop=FALSE]
for (i in 1:ncol(const)) {
cidx <- which(!is.na(const[,i]))
ridx <- match(const[cidx, i], rownames(pars))
pars[cbind(ridx, cidx)] <- colnames(const)[i]
}
- }
+ }
if (shape == "matrix")
out <- pars
if (shape == "vector") {
out <- unique(na.omit(c(t(pars))))
if (ancillary)
- out <- c(out, eqns[fixed])
+ out <- c(out, fixed)
}
out
}
diff --git a/R/model.end.R b/R/model.end.R
index c02e42f..bf0aaba 100644
--- a/R/model.end.R
+++ b/R/model.end.R
@@ -1,4 +1,4 @@
-model.end <- function(res, D) {
+model.end <- function(res, mf) {
res$variance <- -solve(res$hessian)
res$hessian <- NULL
@@ -7,11 +7,11 @@ model.end <- function(res, D) {
res$coefficients <- res$par
res$par <- NULL
- res$terms <- attr(D, "terms")
+ res$terms <- attr(mf, "terms")
- attr(res, "na.message") <- attr(D, "na.message")
- if (!is.null(attr(D, "na.action")))
- res$na.action <- attr(D, "na.action")
+ attr(res, "na.message") <- attr(mf, "na.message")
+ if (!is.null(attr(mf, "na.action")))
+ res$na.action <- attr(mf, "na.action")
res
}
diff --git a/R/model.frame.multiple.R b/R/model.frame.multiple.R
new file mode 100644
index 0000000..44b80a3
--- /dev/null
+++ b/R/model.frame.multiple.R
@@ -0,0 +1,104 @@
+model.frame.multiple <- function (formula,data,...){
+ if(class(formula)[[1]]=="terms"){
+ terms <-formula
+ }else{
+ terms<-terms(formula)
+ }
+ # print(terms)
+ "%w/o%" <- function(x,y) x[!x %in% y]
+ #print("model.frame.multiple is called")
+ toBuildFormula<-function(Xnames,sepp="+"){
+ lng<-length(Xnames)
+ rhs<-NULL
+ if (lng!=0){
+ if(lng==1){
+ rhs=Xnames
+ }else{
+ for (j in 1:(lng-1)){
+ rhs<-paste(rhs,as.name(Xnames[[j]]))
+ rhs<-paste(rhs,sepp)
+ }
+ rhs<-paste(rhs,Xnames[[lng]])
+ }
+ }
+ return (rhs)
+ }
+ eqn<-names(formula)
+ eqn<-attr(terms,"systEqns")
+ nrEquations<-length(eqn)
+ termlabels<-attr(terms,"term.labels")
+ depVars<-attr(terms,"depVars")
+ Xs<-Ys<-tlNew<-dvNew<-list()
+ for (i in 1:nrEquations){
+ rhs<-toBuildFormula(termlabels[[eqn[[i]]]])
+ if(!(is.null(rhs))){
+ rhs<-paste(rhs,"-1")
+ rhs<-as.formula(paste("~",rhs))
+ Xs[[eqn[[i]]]]<-model.matrix.default(rhs,data=data)
+ tlNew[[eqn[[i]]]]<-colnames(Xs[[eqn[[i]]]])
+ tlNew[[eqn[[i]]]]<-gsub("as.factor\\(.*\\)","",tlNew[[eqn[[i]]]],extended=TRUE)
+ colnames(Xs[[eqn[[i]]]])<-tlNew[[eqn[[i]]]]
+ }
+ }
+ depFactors<-attr(terms,"depFactors")
+
+ if(!(is.logical(depFactors)))
+ depVars<- paste("as.factor(",depFactors[[1]],")",sep="")
+ #print(depVars)
+ lhs<-toBuildFormula(unique(unlist(depVars)))
+ if(!(is.null(lhs))){
+ lhs<-paste(lhs,"-1")
+ lhs<-as.formula(paste("~",lhs))
+ Ys<-model.matrix.default(lhs,data=data)
+ dvNew<-colnames(Ys)
+ dvNew<-gsub("as.factor\\(.*\\)","",dvNew,extended=TRUE)
+ colnames(Ys)<-dvNew
+ }
+ attr(terms,"term.labels")[names(tlNew)]<-tlNew
+ attr(terms,"depVars")[names(dvNew)]<-dvNew
+
+ ronames<-rownames(data)
+ ronr<-nrow(data)
+ Xnames<-unique(unlist(tlNew))
+ Ynames<-unique(unlist(dvNew))
+ if(!(is.logical(depFactors)))
+ Ynames<-c(depFactors[[2]],Ynames %w/o% depFactors[[2]])
+ X<-matrix(0,nrow=ronr,ncol=length(Xnames),dimnames=list(ronames,Xnames))
+ Y<-matrix(0,nrow=ronr,ncol=length(Ynames),dimnames=list(ronames,Ynames))
+ if(length(tlNew)>0)
+ for(i in 1:length(tlNew)){
+ xtmp<-intersect(tlNew[[i]],Xnames)
+ X[,xtmp]<-Xs[[i]][,xtmp]
+ }
+ Y<-Ys
+ my.data.frame<-as.data.frame(cbind(Y,X))
+ rhs<-toBuildFormula(Xnames)
+ if(!(is.null(rhs)))
+ rhs<-(paste("~",rhs))
+ else
+ rhs<-"~1"
+ cb<-FALSE
+ if(length(Ynames)>1){
+ lhs<-toBuildFormula(Ynames,",")
+ if (!(is.null(lhs))){
+ lhs<-paste("cbind(",lhs)
+ lhs<-paste(lhs,")")
+ cb<-TRUE
+ }
+ }else{
+ lhs=Ynames
+ }
+ lhs<-as.formula(paste(lhs,rhs))
+#print(lhs)
+#print(names(my.data.frame))
+ Y<-model.frame.default(lhs,data=my.data.frame)
+ result=Y
+ if(cb)
+ names(result)[[1]]<-"response"
+ new.response<-attr(attr(result,"terms"),"response")
+ attr(terms,"response")<-new.response
+ attr(result,"terms")<-terms
+ class(result)<-c(class(result),"multiple")
+ return(result)
+}
+
diff --git a/R/model.matrix.multiple.R b/R/model.matrix.multiple.R
new file mode 100644
index 0000000..28e284e
--- /dev/null
+++ b/R/model.matrix.multiple.R
@@ -0,0 +1,100 @@
+model.matrix.multiple <- function (object,data,shape="compact",eqn=NULL,...){
+
+ intersect <- function(x, y) y[match(x, y, nomatch = 0)]
+
+ toBuildFormula<-function(Xnames,sepp="+"){
+ lng<-length(Xnames)
+ rhs<-NULL
+ if (lng!=0){
+ if(lng==1){
+ rhs=Xnames
+ }else{
+ for (j in 1:(lng-1)){
+ rhs<-paste(rhs,as.name(Xnames[[j]]))
+ rhs<-paste(rhs,sepp)
+ }
+ rhs<-paste(rhs,Xnames[[lng]])
+ }
+ }
+ return (rhs)
+ }
+
+
+ if((shape != "compact") && (shape != "array") && (shape !="stacked"))
+ stop("wrong shape argument! Choose from \"compact\", \"array\" or \"stacked\" \n")
+
+ if(!(any(class(object)=="multiple")))
+ stop("Please run first parse.formula() on you formula ...\n")
+
+ if(!(any(class(data)=="multiple")))
+ data<-model.frame(object,data)
+
+
+ terms<-attr(data,"terms")
+ whiche<-which(eqn %in% names(terms)==FALSE)
+ if (length(whiche)!=0)
+ stop("Unknown eqn name \"",eqn[whiche],"\"\n")
+
+ intercAttr<-attr(terms,"intercept")
+ systEqns<-attr(terms,"systEqns")
+ ancilEqns<-attr(terms,"ancilEqns")
+
+ if (is.null(eqn))
+ eqn=systEqns
+
+ # if (!(all(eqn %in% systEqns)))
+ # stop("all eqn names should be from systematic parameters")
+
+ termlabels<-attr(terms,"term.labels")[eqn]
+ nrEquations<-length(eqn)
+ if (length(eqn)==1)
+ shape="compact"
+
+ Xnames<-unique(unlist(termlabels))
+
+ rhs<-toBuildFormula(Xnames)
+ if(!(is.null(rhs)))
+ rhs<-as.formula(paste("~",rhs))
+ else
+ rhs<-as.formula("~1")
+
+ rawX<-model.matrix.default(rhs,data=data)
+ if (shape=="compact"){
+ result<-rawX
+ if(all(intercAttr==0)){
+ result<-result[,colnames(result)!="(Intercept)"]
+ }
+ attr(terms,"response")<-0
+ attr(result,"terms")<-terms
+ return(result)
+ }
+
+ ronames<-rownames(data)
+ ronr<-nrow(data)
+
+ parsMat<-make.parameters(terms, shape = "matrix", ancillary = FALSE,eqns=eqn)
+ parsVec <- unique(na.omit(c(t(parsMat))))
+
+ result<-list()
+ result<-array(0,dim=c(ronr,length(parsVec),length(eqn)),dimnames=list(ronames,parsVec,eqn))
+ for(i in 1:nrEquations){
+ eqni<-eqn[[i]]
+ whiche<-which(is.na(parsMat[,eqni])==FALSE)
+ result[,,eqni][,parsMat[names(whiche),eqni]]<-rawX[,names(whiche)]
+ }
+
+ if(shape=="array"){
+ res<-result
+ }
+ if(shape=="stacked"){
+ res<-result[,,eqn[[1]]]
+ if(length(eqn)>1)
+ for(i in 2:length(eqn))
+ res<-rbind(res,result[,,eqn[[i]]])
+ rownames(res)<-c(1:nrow(res))
+ }
+ attr(terms,"response")<-0
+ attr(res,"terms")<-terms
+ return(res)
+}
+
diff --git a/R/multi.R b/R/multi.R
new file mode 100644
index 0000000..f92e554
--- /dev/null
+++ b/R/multi.R
@@ -0,0 +1,26 @@
+multi<-function(...){
+ res<-list()
+ mf<-match.call(expand.dots=TRUE)
+ for(i in 2:length(mf)){
+ leveli<-eval(mf[[i]])
+ levelnamei<-names(mf)[[i]]
+ dta<-leveli[[2]]
+ if(class(dta)!="MI")
+ if(is.data.frame(dta[[1]])){
+ for(j in 1:length(dta)){
+ newlevelname<-paste(levelnamei,j,sep="")
+ res[[newlevelname]]<-list(formula=NULL, data=NULL)
+ res[[newlevelname]][[1]]<-leveli[[1]]
+ res[[newlevelname]][[2]]<-leveli[[2]][[j]]
+ }
+ }
+ else{
+ res[[levelnamei]]<-leveli
+ names(res[[levelnamei]])<-c("formula","data")
+ }
+}
+ class(res)<-c("multi", class(res))
+return(res)
+}
+
+
diff --git a/R/parse.formula.R b/R/parse.formula.R
new file mode 100644
index 0000000..2cfc6a9
--- /dev/null
+++ b/R/parse.formula.R
@@ -0,0 +1,323 @@
+parse.formula<-function( formula, model,data=NULL){
+ if(class(formula)[[1]]=="multiple")
+ return(formula)
+ nrUserOpt<-nrUserReq<-nrUserFixed<-nrUserSubreq<-0
+ userOpt<-userReq<-userFixed<-userSubreq<-list()
+
+ fc <- paste("describe.", model, sep = "")
+ if (!exists(fc))
+ stop("describe.",model," does not exsist")
+ modelReq<-do.call(fc,list())
+ modelReq <-modelReq$parameters
+
+ checkNrReq<-function(modelNumEqn,nrUserReq,modelParsReq){
+ if(length(modelNumEqn)==1){
+ if(nrUserReq != modelNumEqn)
+ stop("The parameter \"",modelParsReq,"\" requires ",modelNumEqn, " equation(s).
+ You have provided ", nrUserReq, " See model doc. for more details")
+ }else{
+ if(!(betweenf(nrUserReq,modelNumEqn)))
+ stop("The parameter \"",modelParsReq,"\" requires between ",modelNumEqn[[1]],"
+ and ",modelNumEqn[[2]], " equation(s). You have provided ", nrUserReq, " See model doc. for more details")
+ }
+ }
+
+
+ checkNrOpt<-function(modelNumEqn,nrUserOpt,modelParsOpt,userOpt){
+ if(length(modelNumEqn)==1){
+ if(nrUserOpt > modelNumEqn)
+ stop("The parameter \"",modelParsOpt,"\" requires ",modelNumEqn, " equation(s) with DepVar=FALSE and ExpVar=TRUE. You have provided ", nrUserOpt, " See model doc. for more details")
+ else{
+ if(nrUserOpt<modelNumEqn)
+ for(i in (nrUserOpt+1):modelNumEqn){
+ userOpt[[i]]<-as.formula("~1")
+ if(modelNumEqn==1)
+ names(userOpt)[[i]]<-modelParsOpt
+ else
+ names(userOpt)[[i]]<-paste(modelParsOpt,i,sep="")
+ }
+ }
+ }else{
+ if(!(betweenf(nrUserOpt,modelNumEqn)))
+ if(nrUserOpt < modelNumEqn[[1]])
+ for(i in (nrUserOpt+1):modelNumEqn[[1]]){
+ userOpt[[i]]<-as.formula("~1")
+ names(userOpt)[[i]]<-paste(modelParsOpt,i,sep="")
+ }
+ else
+ stop("The parameter \"",modelParsOpt,"\" requires between ",modelNumEqn[[1]]," and ",modelNumEqn[[2]], " equation(s). You have provided ", nrUserOpt, " See model doc. for more details")
+ }
+ return(userOpt)
+ }
+
+ betweenf<-function(a,range){
+ if (is.finite(range[[2]]))
+ return(a >= range[[1]] && a<=range[[2]])
+ else
+ return(a>=range[[1]])
+ }
+
+ "%w/o%" <- function(x,y) x[!x %in% y]
+
+ matchPars<-function(parName,userNames){
+ res<-c()
+ for(i in 1:length(userNames)){
+ a<-substr(userNames[[i]],nchar(parName)+1,nchar(userNames[[i]]))
+ b<-substr(userNames[[i]],1,nchar(parName))
+ if(b==parName && (!is.na(suppressWarnings(as.numeric(a))) || userNames[[i]]==parName))
+ res<-c(res,userNames[[i]])
+ }
+ return (res)
+ }
+
+ fMode<-function(b){
+ if(b$depVar == TRUE && b$expVar == TRUE) return (1)
+ if(b$depVar == FALSE && b$expVar == TRUE) return (2)
+ if(b$depVar == FALSE && b$expVar == FALSE) return (3)
+ if(b$depVar == TRUE && b$expVar == FALSE) return (4)
+ stop("some error occurred ... please contact the Zelig team")
+ }
+
+ parsType<-lapply(modelReq,fMode)
+ modelParsReq<-names(parsType[parsType==1])
+ modelParsOpt<-names(parsType[parsType==2])
+ modelParsFixed<-names(parsType[parsType==3])
+ modelParsSubreq<-names(parsType[parsType==4])
+
+ modelNrParsReq<-length(modelParsReq)
+ modelNrParsOpt<-length(modelParsOpt)
+ modelNrParsFixed<-length(modelParsFixed)
+ modelNrParsSubreq<-length(modelParsSubreq)
+
+ userNrLevels<-0
+ dataNrLevels<-0
+ userLevels<-c()
+
+ if(class(formula)[[1]]=="formula")
+ formula<-list(formula)
+
+ nreqns <-length(formula)
+
+ if(is.null(names(formula))){
+ if(modelNrParsReq >1)
+ stop("You should name the equations. The model requires more than 1 systematic component. Please see model documentation for more details")
+ for (i in 1:nreqns){
+ eqni<-formula[[i]]
+ if (length(eqni)==3){
+ rootNames<-modelParsReq
+ lhs<-eqni[[2]]
+ rhs<-eqni[[3]]
+ if(length(lhs)>1 && (lhs[[1]]=="cbind" || lhs[[1]]=="as.factor" || lhs[[1]]=="id")){
+ if( lhs[[1]]=="cbind"){
+ rhs=deparse(rhs)
+ g<- as.list(lhs)[-1]
+ for (j in 1:length(g)){
+ e<-paste(g[[j]],"~",sep="")
+ if(rhs!="1"){
+ nrUserReq=nrUserReq+1
+ userReq[[nrUserReq]]<-as.formula(paste(e,rhs,sep=""))
+ }else{
+ nrUserSubreq=nrUserSubreq+1
+ userSubreq[[nrUserSubreq]]<-as.formula(paste(e,rhs,sep=""))
+ }
+ }
+ }else{
+ if(is.null(data))
+ stop("Data argument is required when you use as.factor() or id() as a dependent variable\n")
+ if(lhs[[1]]=="as.factor"){
+ varname<-as.character(lhs[[2]])
+ userLevels<-levels(as.factor(data[[varname]]))[-1]
+ userNrLevels<-length(userLevels)
+ for (j in 1:userNrLevels){
+ e<-paste("id(",lhs[[2]],",\"",userLevels[[j]],"\")","~",sep="")
+ if(rhs!="1"){
+ nrUserReq=nrUserReq+1
+ userReq[[nrUserReq]]<-as.formula(paste(e,deparse(rhs),sep=""))
+ }else{
+ nrUserSubreq=nrUserSubreq+1
+ userSubreq[[nrUserSubreq]]<-as.formula(paste(e,deparse(rhs),sep=""))
+ }
+ }
+ }else{
+ varname<-as.character(lhs[[2]])
+ userLevels<-c(userLevels,lhs[[3]])
+ userNrLevels<-length(userLevels)
+ levels<-levels(data[[varname]])
+ lhs<-deparse(lhs)
+ rhs<-deparse(rhs)
+ e<-paste(lhs,"~",sep="")
+ if(rhs !="1"){
+ nrUserReq=nrUserReq+1
+ userReq[[nrUserReq]]<-as.formula(paste(e,rhs,sep=""))
+ }else{
+ nrUserSubreq<-nrUserSubreq+1
+ userSubreq[[nrUserSubreq]]<-as.formula(paste(e,rhs,sep=""))
+ }
+ }
+ }
+ }else{
+ lhs<-deparse(lhs)
+ rhs<-deparse(rhs)
+ e<-paste(lhs,"~",sep="")
+ if(rhs !="1"){
+ nrUserReq=nrUserReq+1
+ userReq[[nrUserReq]]<-as.formula(paste(e,rhs,sep=""))
+ }else{
+ nrUserSubreq<-nrUserSubreq+1
+ userSubreq[[nrUserSubreq]]<-as.formula(paste(e,rhs,sep=""))
+ }
+ }
+ }else{
+ rhs<-deparse(eqni[[2]])
+ if(rhs !="1"){
+ nrUserOpt=nrUserOpt+1
+ userOpt[[nrUserOpt]]<-as.formula(paste("~",rhs,sep=""))
+ }else{
+ nrUserFixed=nrUserFixed+1
+ userFixed[[nrUserFixed]]<-as.formula(paste("~",rhs,sep=""))
+ }
+ }
+ }
+ if (modelNrParsOpt==0){
+ if (nrUserOpt !=0){
+ stop("the equation(s) ",userOpt," does not match model requirements!")}
+ }else{
+ modelNumEqn<-modelReq[[modelParsOpt]]$equations
+ userOpt<-checkNrOpt(modelNumEqn,nrUserOpt,modelParsOpt,userOpt)
+ if(length(userOpt)==1)
+ names(userOpt)<-modelParsOpt
+ else
+ names(userOpt)<-paste(modelParsOpt,1:length(userOpt),sep="")
+ }
+
+ if(length(modelParsFixed)>0){
+ modelNumFixedEqns<-modelReq[[modelParsFixed]]$equations
+ for(i in 1:modelNumFixedEqns)
+ userFixed[[i]]<-as.formula("~1")
+ if(modelNumFixedEqns==1)
+ names(userFixed)<-modelParsFixed
+ else
+ names(userFixed)<-paste(modelParsFixed,1:modelNumFixedEqns,sep="")
+ }
+ if (modelNrParsReq==0){
+ if (nrUserReq !=0){
+ stop("the equation(s) ",userReq," does not match model requirements!")}
+ }else{
+ modelNumEqn<-modelReq[[modelParsReq]]$equations
+ checkNrReq(modelNumEqn,nrUserReq,modelParsReq)
+ if(userNrLevels>0){
+ if(userNrLevels !=nrUserReq)
+ stop("The number of equation for the systematic component should be equal to the number of levels -1\n")
+ names(userReq)<-userLevels
+ }else{
+ if(nrUserReq==1)
+ names(userReq)<-modelParsReq
+ else
+ names(userReq)<-paste(modelParsReq,1:length(userReq),sep="")
+ }
+ }
+
+ if (modelNrParsSubreq==0){
+ if (nrUserSubreq !=0){
+ stop("the equation(s) ",userSubreq," does not match model requirements!")}
+ }else{
+ modelNumEqn<-modelReq[[modelParsSubreq]]$equations
+ checkNrReq(modelNumEqn,nrUserSubreq,modelParsSubreq)
+ if(nrUserSubreq==1)
+ names(userSubreq)<-modelParsSubreq
+ else
+ names(userSubreq)<-paste(modelParsSubreq,1:length(userSubreq),sep="")
+ }
+ result<-c(userReq,userOpt,userFixed,userSubreq)
+ }else{
+ modelPars<-names(modelReq)
+ parsS<-names(sort(sapply(modelPars,nchar),decreasing=TRUE))
+ userNames<-names(formula)
+ userEqnNamesByPars<-list()
+ tmpUserNames<-userNames
+ for (i in 1:length(parsS)){
+ userEqnNamesByPars[[parsS[[i]]]]<-matchPars(parsS[[i]],tmpUserNames)
+ tmpUserNames<-"%w/o%"(tmpUserNames,userEqnNamesByPars[[parsS[[i]]]])
+ }
+ tmp<-"%w/o%"(userNames,unlist(userEqnNamesByPars))
+ if (length(tmp)>0)
+ stop("Ambigous equation name ","\"",tmp,"\"")
+ res<-list()
+ userPars<-names(userEqnNamesByPars)
+ for (i in 1:length(modelPars)){
+ modelPar<-modelPars[[i]]
+ userNumEqn<-length(userEqnNamesByPars[[modelPar]])
+ modelNumEqn<-modelReq[[modelPar]]$equations
+ mode<-fMode(modelReq[[modelPar]])
+ tmplst<-formula[userEqnNamesByPars[[modelPar]]]
+ if(length(modelNumEqn)==1 && modelNumEqn==1)
+ tmpNames<-modelPar
+ else
+ tmpNames<-paste(modelPar,1:userNumEqn,sep="")
+ if(mode==1){
+ whiche<-which(lapply(formula[(userEqnNamesByPars[[modelPar]])],length)!=3)
+ if(length(whiche)!=0)
+ stop("The equation ",formula[[names(whiche)]]," is not conform model requirements or its name is ambigous . DepVar/ExpVar is missing.\n")
+ checkNrReq(modelNumEqn,userNumEqn,modelPar)
+ whiche<-which((names(tmplst) %in% tmpNames)==FALSE)
+ if(length(whiche)!=0){
+ warning("The name \"",names(tmplst)[whiche],"\" is ambigous. The equations of the paramter \"",modelPar,"\" are renamed\n")
+ names(tmplst)<-tmpNames
+ }
+ }else{
+ if(mode==2){
+ whiche<-which(lapply(formula[(userEqnNamesByPars[[modelPar]])],length)!=2)
+ if(length(whiche)!=0)
+ stop("The equation ",formula[names(whiche)]," is not conform model requirements or its name is ambigous A .\n")
+ whiche<-which((names(tmplst) %in% tmpNames)==FALSE)
+ if(length(whiche)!=0){
+ warning("The name \"",names(tmplst)[whiche],"\" is ambigous. The equations of the paramter \"",modelPar,"\" are renamed\n")
+ names(tmplst)<-tmpNames
+ }
+ tmplst<- checkNrOpt(modelNumEqn,userNumEqn,modelPar,tmplst)
+ }else{
+ if (mode==3){
+ whiche<-which(tmplst !="~1")
+ if(length(whiche)>0)
+ warning("You cannot specify a formula for the parameter \"",modelPar,"\" . All your equation for this parameter are set to their default value.For example your equation:\n",deparse(formula[names(whiche)]),"\n")
+ if(userNumEqn !=modelNumEqn)
+ warning("The parameter \"",modelPar,"\" requires ",modelNumEqn, "equation(s). You are providing ",userNumEqn, " equation(s) for this parameter. This problem is fixed. All the equations for this parameter are set to the default value \n")
+
+ tmplst<-list()
+ if(modelNumEqn==1)
+ tmpname<-modelPar
+ else
+ tmpname<-paste(modelPar,1:modelNumEqn,sep="")
+ for(i in 1:modelNumEqn)
+ tmplst[[tmpname[[i]]]]<-as.formula("~1")
+ }else{
+ if(mode==4)
+ {
+ whiche<-which(lapply(formula[(userEqnNamesByPars[[modelPar]])],length)!=3)
+ whicha<-which(lapply(formula[(userEqnNamesByPars[[modelPar]])],FUN=function(a){if (a[[3]]=="1") return (TRUE) else return(FALSE)})==FALSE)
+ if(length(whiche)!=0 )
+ stop("The equation ",formula[names(whiche)]," is not conform model requirements or its name is ambigous . DepVar/ExpVar is missing.\n")
+ else{
+ if (length(whicha)!=0)
+ stop("The equation ",formula[names(whicha)]," is not conform model requirements or its name is ambigous . Its right hand side shoule be \"1\".\n")
+ }
+ checkNrReq(modelNumEqn, userNumEqn, modelPar)
+ whiche<-which((names(tmplst) %in% tmpNames)==FALSE)
+ if(length(whiche)!=0){
+ warning("The name \"",names(tmplst)[whiche],"\" is ambigous. The equations of the paramter \"",modelPar,"\" are renamed\n")
+ names(tmplst)<-tmpNames
+ }
+ }
+ }
+ }
+ }
+ res[[modelPar]]<-tmplst
+ }
+ result<-c()
+ for(i in 1:length(res))
+ result<-c(result,res[[i]])
+
+ }
+ class(result)<-c("multiple","list")
+ return(result)
+}
diff --git a/R/parse.par.R b/R/parse.par.R
new file mode 100644
index 0000000..65abc39
--- /dev/null
+++ b/R/parse.par.R
@@ -0,0 +1,50 @@
+parse.par <- function(par, terms,shape = "matrix", eqn=NULL) {
+ "%w/o%" <- function(x,y) x[!x %in% y]
+ if (is.null(shape)) {
+ if (any(class(terms) == "multiple"))
+ shape <- "matrix"
+ else
+ shape <- "vector"
+ }
+ if(is.null(eqn))
+ eqn<-attr(terms,"systEqns")
+ if (!shape %in% c("matrix", "vector"))
+ stop("not a valid 'shape' for parameters. Choose from \"matrix\" or \"vector\".")
+ if (any(class(terms) == "multiple")) {
+ allidx <- make.parameters(terms = terms, shape = "vector")
+ idx <- make.parameters(terms = terms,eqn=eqn, shape = "vector")
+ mat <- t(make.parameters(terms = terms,eqn=eqn, shape = "matrix"))
+ if(length(par)==length(allidx))
+ par.names<-allidx
+ else
+ par.names<-idx
+ ancil <-attr(terms,"ancilEqns")
+ syst<-attr(terms,"systEqns")
+ if (length(syst) == 1)
+ shape <- "vector"
+ if (any(eqn %in% ancil)) {
+ if (any(eqn %in% syst)) {
+ stop(" eqn cannot include both systematic and ancillary \n parameters at the same time.")
+ }
+ else
+ ret.val <- par[par.names %in% idx]
+ }
+ else { ## if eqn to be returned is a systematic component
+ subs<-mat
+ out <- matrix(0, nrow = nrow(subs), ncol = ncol(subs),
+ dimnames = dimnames(subs))
+ for(i in 1:nrow(out))
+ for(j in 1:ncol(out))
+ if(!is.na(subs[i,j]))
+ out[i,j]<-par[par.names %in% subs[i,j]]
+ if (shape == "matrix")
+ ret.val <- t(out)
+ else {
+ ret.val <- par[par.names %in% idx]
+ }
+ }
+ }
+ ret.val
+}
+
+
diff --git a/R/MULTIPLE_NEW/put.start.R b/R/put.start.R
similarity index 89%
rename from R/MULTIPLE_NEW/put.start.R
rename to R/put.start.R
index d59ac03..fd04392 100644
--- a/R/MULTIPLE_NEW/put.start.R
+++ b/R/put.start.R
@@ -1,9 +1,8 @@
put.start <- function(start.val, value, terms, eqn) {
- if (!any(class(terms) == "list"))
+ if (!any(class(terms) == "multiple"))
stop("'put.start()' works with 'parse.formula()'. Use that first!")
idx <- names(start.val)
const <- attr(terms, "constraints")
- const[const == "0"] <- NA
if (!is.logical(const)) {
for (var in colnames(const)) {
eqns <- paste(names(na.omit(const[,var])), collapse = ":")
diff --git a/R/qi.MCMCZelig.R b/R/qi.MCMCZelig.R
index c429046..4b887d0 100644
--- a/R/qi.MCMCZelig.R
+++ b/R/qi.MCMCZelig.R
@@ -33,7 +33,7 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
else if (model =="normal.bayes") {
nvar <- ncol(object$coefficients)
coef <- object$coefficients[,1:(nvar-1)]
- qi$ev <- coef %*% t(x)
+ qi$ev <- ev <- coef %*% t(x)
qi$pr <- rnorm(nrow(qi$ev), qi$ev,
sqrt(object$coefficients[,nvar]))
qi.name <- list(ev = "Expected Values: E(Y|X)", pr = "Predicted Values:Y|X")
@@ -49,11 +49,13 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
L2 <- (object$above-eta)/sig
L1 <- (object$below-eta)/sig
##cev <- eta + sig*(dnorm(L1)-dnorm(L2))/(pnorm(L2)-pnorm(L1))
+ temp1 <- pnorm(L1)*object$below
if (object$below==-Inf) temp1<-0
- else temp1 <- pnorm(L1)*object$below
+
+ temp2 <- (1-pnorm(L2))*object$above
if (object$above==Inf) temp2<-0
- else temp2 <- (1-pnorm(L2))*object$above
- qi$ev <- temp1+eta*(pnorm(L2)-pnorm(L1))+sig*(dnorm(L1)-dnorm(L2))+temp2
+
+ qi$ev <-ev <- temp1+eta*(pnorm(L2)-pnorm(L1))+sig*(dnorm(L1)-dnorm(L2))+temp2
qi.name <- list(ev = "Expected Values: E(Y|X)")
}
else if (model == "poisson.bayes") {
@@ -87,17 +89,19 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
}
else if (model == "normal.bayes") {
ev1 <- eta1
- qi$fd <-ev1-ev
+ qi$fd <- ev1 - ev
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))
+ temp1 <- pnorm(L1)*object$below
if (object$below==-Inf) temp1<-0
- else temp1 <- pnorm(L1)*object$below
+
+ temp2 <- (1-pnorm(L2))*object$above
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
qi$fd <-ev1-ev
qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
diff --git a/R/qi.survreg.R b/R/qi.survreg.R
index f22259d..26a9b1e 100644
--- a/R/qi.survreg.R
+++ b/R/qi.survreg.R
@@ -8,13 +8,16 @@ qi.survreg <- function(object, simpar, x, x1 = NULL, y = NULL) {
else
sim.scale <- rep(object$scale, nrow(simpar))
}
- else if (model == "lognorm")
+ else if (model %in% c("lognorm", "tobit"))
sim.scale <- simpar[,(k+1):ncol(simpar)]
if (!is.null(y)) {
status <- y[,2]
y <- y[,1]
}
- link <- survreg.distributions[[object$dist]]$itrans
+ if (model %in% c("weibull", "Weibull", "lognorm", "exp"))
+ link <- survreg.distributions[[object$dist]]$itrans
+ else if (model == "tobit")
+ link <- function(x) x
ev.surv <- function(model, sim.coef, sim.scale, x, link) {
eta <- sim.coef %*% t(x)
theta <- as.matrix(apply(eta, 2, link))
@@ -22,11 +25,11 @@ qi.survreg <- function(object, simpar, x, x1 = NULL, y = NULL) {
ev <- exp(log(theta) + 0.5*(exp(sim.scale))^2)
dimnames(ev) <- dimnames(theta)
}
- else if (model == "weibull" || model == "Weibull") {
+ else if (model %in% c("weibull", "Weibull")) {
ev <- theta * gamma(1 + exp(sim.scale))
dimnames(ev) <- dimnames(theta)
}
- else if (model == "exp") {
+ else if (model %in% c("exp", "tobit")) {
ev <- theta
}
list(ev = as.matrix(ev), theta = as.matrix(theta))
@@ -34,12 +37,15 @@ qi.survreg <- function(object, simpar, x, x1 = NULL, y = NULL) {
pr.surv <- function(model, theta, sim.scale, ev) {
if (model == "exp")
pr <- rexp(length(ev), rate = 1/ev)
- else if (model == "weibull" || model == "Weibull")
+ else if (model %in% c("weibull", "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))
+ else if (model == "tobit") {
+ pr <- rnorm(length(ev), mean = ev, sd = exp(sim.scale))
+ }
pr
}
ev <- ev.surv(model, sim.coef, sim.scale, x, link)
diff --git a/R/relogit.R b/R/relogit.R
index 6729282..6052c42 100644
--- a/R/relogit.R
+++ b/R/relogit.R
@@ -15,6 +15,8 @@ relogit <- function(formula, data = sys.parent(), tau = NULL,
else
weighting <- TRUE
}
+ else
+ weighting <- FALSE
if (length(tau) > 2)
stop("tau must be a vector of length less than or equal to 2")
else if (length(tau)==2) {
@@ -79,9 +81,13 @@ relogit <- function(formula, data = sys.parent(), tau = NULL,
res$coefficients["(Intercept)"] <- res$coefficients["(Intercept)"] -
log(((1-tau)/tau) * (ybar/(1-ybar)))
res$prior.correct <- TRUE
+ res$weighting <- FALSE
}
else
res$prior.correct <- FALSE
+ if (is.null(res$weighting))
+ res$weighting <- FALSE
+
res$linear.predictors <- t(res$coefficients) %*% t(X)
res$fitted.values <- 1/(1+exp(-res$linear.predictors))
res$zelig <- "relogit"
diff --git a/R/MULTIPLE_NEW/set.start.R b/R/set.start.R
similarity index 82%
rename from R/MULTIPLE_NEW/set.start.R
rename to R/set.start.R
index 18b370b..ff4c721 100644
--- a/R/MULTIPLE_NEW/set.start.R
+++ b/R/set.start.R
@@ -1,5 +1,5 @@
-set.start <- function(start.val, terms) {
- if (any(class(terms) == "list"))
+set.start <- function(start.val = NULL, terms) {
+ if (any(class(terms) == "multiple"))
labs <- make.parameters(terms = terms, shape = "vector", ancillary = TRUE)
else
labs <- attr(terms, "term.labels")
diff --git a/R/setx.default.R b/R/setx.default.R
index a0ddbda..94d340d 100644
--- a/R/setx.default.R
+++ b/R/setx.default.R
@@ -1,7 +1,7 @@
setx.default <- function(object, fn = list(numeric = mean, ordered =
median, other = mode), data = NULL,
cond = FALSE, counter = NULL, ...){
- mf <- match.call()
+ mc <- match.call()
mode <- function(x){
tb <- tapply(x, x, length)
if(is.factor(x))
@@ -50,19 +50,30 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
stop("min cannot be calculated for this data type")
return(value)
}
- if (any(class(object)=="vglm"))
- tt <- object at terms$terms
+ tt <- terms(object)
+ tt.attr <- attributes(tt)
+ env <- tt.attr$.Environment
+ if (is.null(env))
+ env <- parent.frame()
+ ## original data
+ if (is.null(data))
+ if (is.data.frame(object$data))
+ dta <- object$data
+ else
+ dta <- eval(object$call$data, envir = env)
else
- tt <- terms(object)
- if (is.null(data))
- odta <- eval(object$call$data, sys.parent())
+ dta <- as.data.frame(data)
+ ## extract variables we need
+ mf <- model.frame(tt, data = dta, na.action = na.pass)
+ vars <- all.vars(object$call)
+ if (!is.null(tt.attr$response) && tt.attr$response)
+ resvars <- all.vars(tt.attr$variables[[1+tt.attr$response]])
else
- odta <- data
- data <- dta <- odta[attributes(model.frame(tt, odta))$row.names,]
- vars <- names(dta)
+ resvars <- NULL
+ data <- dta[complete.cases(mf), names(dta)%in%vars]
if (!is.null(counter)) {
if (!any(counter == vars))
- stop(paste("the variable specified for counter is not used in the model"))
+ stop("the variable specified for counter is not used in the model")
treat <- data[, names(data)==counter]
if(is.numeric(treat)) {
data[treat==1, names(data)==counter] <- 0
@@ -76,7 +87,7 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
data[treat==0, names(data)==counter] <- lev[2]
}
else
- stop(paste("counter only takes a binary variable"))
+ stop("counter only takes a binary variable")
}
else if(is.logical(treat)) {
treat <- as.numeric(treat)
@@ -84,15 +95,17 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
data[treat==0, names(data)==counter] <- TRUE
}
else
- stop(paste("not supported variable type for counter"))
+ stop("not supported variable type for counter")
if(!cond)
- stop(paste("if counter is specified, cond must be TRUE"))
+ stop("if counter is specified, cond must be TRUE")
}
if (cond) {
if (is.null(data))
- stop(paste("if cond = TRUE, you must specify the data frame."))
+ stop("if cond = TRUE, you must specify the data frame.")
+ if (is.null(mc$fn))
+ fn <- NULL
if (!is.null(fn)) {
- warning(paste("when cond = TRUE, fn is coerced to NULL"))
+ warning("when cond = TRUE, fn is coerced to NULL")
fn <- NULL
}
}
@@ -117,19 +130,21 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
fn$other <- mode
}
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
+ if (!(colnames(data)[i] %in% resvars)) {
+ 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))]
+ opt <- vars[na.omit(pmatch(names(mc), vars))]
maxl <- 1
if (length(opt) > 0)
for (i in 1:length(opt)) {
- value <- eval(mf[[opt[i]]], sys.parent())
+ value <- eval(mc[[opt[i]]], envir = env)
lv <- length(value)
if (lv>1)
if (maxl==1 || maxl==lv) {
@@ -154,7 +169,7 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
names(data) <- vars
}
if (cond) {
- X <- model.frame(tt, odta)
+ X <- model.frame(tt, data = dta)
if (!is.null(counter)) {
X <- list(treat=X[treat==1,], control=X[treat==0,])
class(X$treat) <- class(X$control) <- c("data.frame", "cond")
@@ -164,7 +179,7 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
class(X) <- c("data.frame", "cond")
}
else
- X <- as.data.frame(model.matrix(delete.response(tt), data = data))
+ X <- as.data.frame(model.matrix(tt, data = data))
return(X)
}
diff --git a/R/sim.default.R b/R/sim.default.R
index f81c0b6..a0cc2c3 100644
--- a/R/sim.default.R
+++ b/R/sim.default.R
@@ -36,6 +36,7 @@ sim.default <- function(object, x, x1=NULL, num=c(1000, 100),
}
simqi <- qi(object, simpar = simpar, x = x, x1 = x1, y = NULL)
c <- match.call()
+ c[[1]] <- as.name("sim")
c$num <- num
res <- list(x=x, x1=x1, call = c, zelig.call = object$call,
par = simpar, qi=simqi$qi, qi.name=simqi$qi.name)
diff --git a/R/sim.setx.MI.R b/R/sim.setx.MI.R
index 11144b0..ae3ccbe 100644
--- a/R/sim.setx.MI.R
+++ b/R/sim.setx.MI.R
@@ -11,7 +11,8 @@ sim.setx.MI <- function(object, x, x1 = NULL, num = c(1000, 100), prev = NULL,
ca <- match.call()
if (!any(class(x) == "cond")) {
simpar <- MIsimulation(object, num, prev, bootstrap, bootfn, ...)
- simqi <- qi(object[[1]], simpar = simpar, x = x, x1 = x1)
+ simqi <- qi(object[[1]], simpar = simpar, x = as.matrix(x),
+ x1 = if (!is.null(x1)) as.matrix(x1))
ca$num <- num
res <- list(x = x, x1 = x1, call = ca, zelig.call = object[[1]]$call,
par = simpar, qi = simqi$qi, qi.name = simqi$qi.name)
diff --git a/R/summary.glm.robust.R b/R/summary.glm.robust.R
index 8556a98..ea731fe 100644
--- a/R/summary.glm.robust.R
+++ b/R/summary.glm.robust.R
@@ -1,17 +1,16 @@
summary.glm.robust <- function(object, ...) {
- class(object) <- "glm"
+ class(object) <- c("glm", "lm")
res <- summary.glm(object, ...)
if (is.null(object$robust)) {
res$cov.unscaled <- covmat.unscaled <- vcovHAC(object)
res$robust <- "vcovHAC"
- }
- else {
+ } 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 <- covmat.unscaled <- eval(do.call(fn, arg))
+ res$cov.unscaled <- covmat.unscaled <- eval(do.call(fn, args=arg))
}
res$cov.scaled <- covmat <- covmat.unscaled*res$dispersion
if (!is.null(res$correlation)) {
diff --git a/R/systemfit/2sls.tex b/R/systemfit/2sls.tex
new file mode 100644
index 0000000..890ad21
--- /dev/null
+++ b/R/systemfit/2sls.tex
@@ -0,0 +1 @@
+\documentclass[12pt]{book}%
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{fullpage}
\usepackage{multirow}
\usepackage{hyperref}%
\usepackage{bibentry}
\begin{document}
\subsection{\texttt{2sls}: Two Stage Least Squares}
\label{2sls}
\texttt{2sls} provides consistent estimates for linear regression models with
some explanatory variable (the instrumental variable)
correlated with the error term.
In this situation, ordinary least squares fails to provide consistent
estimates. The name two-stage least squares stems from the two regressions
in the estimation procedure. In stage one, an ordinary least squares
prediction of the instrumental variable is obtained from regressing it on
the instrument variables. In stage two, the coefficients of interst are
estimated using ordinary least square after substituting the instrumental
variable by its predictions from stage one.
\subsubsection{Syntax}
\begin{verbatim}
> fml <- list ("mu" = Y ~ X + Z,
"inst" = Z ~ W + X)
> z.out <- zelig(formula = fml, model = "2sls", data = mydata)
> x.out <- setx(z.out)
> s.out <- sim(z.out, x = x.out)
\end{verbatim}
\subsubsection{Inputs}
\texttt{2sls} regression take the following inputs:
\begin{itemize}
\item \texttt{formula}:a list of the main equation and instrumental variable
equation. The first object in the list \texttt{mu} corresponds to the
regression model needs to be estimated. The second list object \texttt{inst}
specifies the regression model for the instrumental variable \texttt{Z}.
For example:
\begin{verbatim}
> fml <- list ("mu" = Y ~ X + Z,
"inst" = Z ~ W + X)
\end{verbatim}
\begin{itemize}
\item \texttt{Y}: the dependent variable of interest.
\item \texttt{Z}: the instrumental variable.
\item \texttt{W}: exogenous instrument variables.
\end{itemize}
\end{itemize}
\subsubsection{Additional Inputs}
\texttt{2sls} takes the following additional inputs for model
specifications:
\begin{itemize}
\item \texttt{TX}: an optional matrix to transform the regressor
matrix and, hence, also the coefficient vector (see \ref{details}). Default is \texttt{NULL}.
\item \texttt{rcovformula}: formula to calculate the estimated residual covariance
matrix (see \ref{details}). Default is equal to 1.
\item \texttt{probdfsys}: use the degrees of freedom of the whole system
(in place of the degrees of freedom of the single equation to calculate probability
values for the t-test of individual parameters.
\item \texttt{single.eq.sigma}: use different $\sigma^2$ for each single
equation to calculate the covariance matrix and the standard errors of the coefficients.
\item \texttt{solvetol}: tolerance level for detecting linear dependencies when
inverting a matrix or calculating a determinant. Default is \texttt {solvetol}=.Machine\$double.eps.
\item \texttt{saveMemory}: logical. Save memory by omitting some calculation that are
not crucial for the basic estimate (e.g McElroy's $R^2$).
\end{itemize}
\subsubsection{Details}
\label{details}
\begin{itemize}
\item \texttt{TX}: The matrix \texttt{TX} transforms the regressor matrix
($X$) by $X\ast=X \times TX$. Thus,
the vector of coefficients is now $b=TX \times b\ast$ where $b$ is the
original(stacked)
vector of all coefficients and $b\ast$ is the new coefficient vector
that is estimated instead.
Thus, the elements of vector $b$ and $b_i = \sum_j TX_{ij}\times b_j\ast$. The $TX$ matrix can be
used to change the order of the coefficients and also to restrict coefficients (if $TX$ has
less columns than it has rows).
\item \texttt{rcovformula}: The formula to calculate the estimated covariance matrix of the residuals($\hat{\Sigma}$)can be one
of the following (see Judge et al., 1955, p.469):
if \texttt{rcovformula}= 0:
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T}
\end{eqnarray*}
if \texttt{rcovformula}= 1 or \texttt{rcovformula}='geomean':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{\sqrt{(T-k_i)\times (T-k_j)}}
\end{eqnarray*}
if \texttt{rcovformula}= 2 or \texttt{rcovformula}='Theil':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T-k_i-k_j+tr[X_i(X_i\prime X_i)^{-1}X_i\prime X_j(X_j\prime X_j)^{-1}X_j\prime]}
\end{eqnarray*}
if \texttt{rcovformula}= 3 or \texttt{rcovformula}='max':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T-max(k_i,k_j)}
\end{eqnarray*}
If $i = j$, formula 1, 2, and 3 are equal. All these three formulas yield unbiased estimators
for the diagonal elements of the residual covariance matrix. If $i neq j$, only formula 2
yields an unbiased estimator for the residual covariance matrix, but it is not necessarily
positive semidefinit. Thus, it is doubtful whether formula 2 is really superior to formula 1
\end{itemize}
\subsubsection{Examples}
\subsubsection{Model}
Let's consider the following regression model,
\begin{eqnarray*}
Y_i=X_i\beta + Z_i\gamma + \epsilon_i, \quad i=1,\ldots,N
\end{eqnarray*}
where $Y_i$ is the dependent variable,
$X_i = (X_{1i},\ldots, X_{Ni})$ is the vector of explanatory variables,
$\beta$ is the vector of coefficients of the explanatory variables $X_i$,
$Z_i$
is the problematic explanatory variable, and $\gamma$ is the coefficient
of $Z_i$. In the equation, there is a direct dependence of $Z_i$
on the structural disturbances of $\epsilon$.
\begin{itemize}
\item The \emph{stochastic component} is given by
\begin{eqnarray*}
\epsilon_i & \sim & {\cal N}(0, \sigma^2), \quad {\rm and} \quad
{\rm cov}(Z_i, \epsilon_i) \ne 0,
\end{eqnarray*}
\item The \emph{systematic component} is given by:
\begin{eqnarray*}
\mu_{i}= E(Y_i)= X_{i}\beta + Z_i\gamma,
\end{eqnarray*}
\end{itemize}
\noindent To correct the problem caused by the correlation of $Z_i$ and $\epsilon$,
two stage least squares utilizes two steps:
\begin{itemize}
\item \emph{Stage 1}: A new instrumental variable $\hat{Z}$ is created
for $Z_i$ which is the ordinary least squares predictions from regressing
$Z_i$ on a set of exogenous instruments $W$ and $X$.
\begin{eqnarray*}
\widehat{Z_i} = \widetilde{W}_i[(\widetilde{W}^\top\widetilde{W})^{-1}\widetilde{W}^\top Z]
\end{eqnarray*}
where $\widetilde{W} = (W,X)$
\item \emph{Stage 2}: Substitute for $\hat{Z}_i$ for $Z_i$ in the original
equation, estimate $\beta$ and $\gamma$ by ordinary least squares regression
of $Y$ on $X$ and $\hat{Z}$ as in the following equation.
\begin{eqnarray*}
Y_i=X_i\beta + \widehat{Z_i}\gamma + \epsilon_i, \quad {\rm for}
\quad i=1,\ldots,N
\end{eqnarray*}
\end{itemize}
\subsubsection{See Also}
For information about three stage least square regression, see
\Sref{3sls} and \texttt{help(3sls)}.
For information about seemingly unrelated regression, see
\Sref{sur} and \texttt{help(sur)}.
\subsubsection{Quantities of Interest}
\subsubsection{Output Values}
The output of each Zelig command contains useful information which you may
view. For example, if you run:
\begin{verbatim}
z.out <- zelig(formula=fml, model = "2sls", data)
\end{verbatim}
\noindent then you may examine the available information in \texttt{z.out} by
using \texttt{names(z.out)}, see the draws from the posterior distribution of
the \texttt{coefficients} by using \texttt{z.out\$coefficients}, and view a default
summary of information through \texttt{summary(z.out)}. Other elements
available through the \texttt{\$} operator are listed below:
\begin{itemize}
\item \texttt{h}: matrix of all (diagonally stacked) instrumental variables.
\item \texttt{single.eq.sigma}: different $\sigma^2$s for each single equation?.
\end{itemize}
\input{systemfit_output}
\begin{itemize}
\item \texttt{inst*}: instruments of the ith equation.
\item \texttt{h*}: matrix of instrumental variables of the ith equation.
\end{itemize}
\subsubsection{Contributors}
The \texttt{2sls} function is adapted from \input{contributors_systemfit}
\end{document}
\ No newline at end of file
diff --git a/R/systemfit/3sls.tex b/R/systemfit/3sls.tex
new file mode 100644
index 0000000..bb266b3
--- /dev/null
+++ b/R/systemfit/3sls.tex
@@ -0,0 +1 @@
+\documentclass[12pt]{book}%
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{fullpage}
\usepackage{multirow}
\usepackage{hyperref}%
\usepackage{bibentry}
\begin{document}
\subsection{\texttt{3sls}: Three Stage Least Squares}
\label{3sls}
\texttt{3sls} is a combination of two stage least squares
and seemingly unrelated regression. It provides consistent estimates for linear regression models with
explanatory variables correlated with the error term. It also extends ordinary least squares
analysis to estimate system of linear equations with correlated error terms
\subsubsection{Syntax}
\begin{verbatim}
> fml <- list ("mu1" = Y1 ~ X1 + Z1,
"mu2" = Y2 ~ X2 + Z2,
"inst1" = Z1 ~ W1 + X1,
"inst2" = Z2 ~ W2 + X2)
> z.out <- zelig(formula = fml, model = "3sls", data = mydata)
> x.out <- setx(z.out)
> s.out <- sim(z.out, x = x.out)
\end{verbatim}
\subsubsection{Inputs}
\texttt{3sls} regression specification requires at least two sets of equations. The first set of $M$ euqations
corresponds to the $M$ dependent variables ($Y_1,\ldots,Y_M$) to be estimated. The second set of equations ($Z$)
corresponds to the instrumental variables in the $M$ equations.
\begin{itemize}
\item \texttt{formula}:a list of the system of equations and instrumental variable
equations. The system of equations is listed first as \texttt{mu}s. The equations
for the instrumental variables are listed next as \texttt{inst}s.
For example:
\begin{verbatim}
> fml <- list ("mu1" = Y1 ~ X1 + Z1,
"mu2" = Y2 ~ X2 + Z2,
"inst1" = Z1 ~ W1 + X1,
"inst2" = Z2 ~ W2 + X2)
\end{verbatim}
\texttt{"mu1"} is the first equation in the two equation model with \texttt{Y1}
as the dependent variable and \texttt{X1} and \texttt{Z1} as
the explanatory variables. \texttt{"mu2"} is the second equation with
\texttt{Y2} as the dependent variable
and \texttt{X2} and \texttt{Z2} as the explanatory variables.
\texttt{Z1} and \texttt{Z2} are also problematic endogenous variables, so
they are estimated through instruments in the \texttt{"inst1"}
and \texttt{"inst2"} equations.
\item \texttt{Y}: dependent variables of interest in the system of equations.
\item \texttt{Z}: the problematic explanatory variables correlated with
the error term.
\item \texttt{W}: exogenous instrument variables used to estimate the
problematic explanatory variables (\texttt{Z})
\end{itemize}
\subsubsection{Additional Inputs}
\texttt{3sls} takes the following additional inputs for model
specifications:
\begin{itemize}
\item \texttt{TX}: an optional matrix to transform the regressor
matrix and, hence, also the coefficient vector (see details). Default is \texttt{NULL}.
\item \texttt{maxiter}: maximum number of iterations.
\item \texttt{tol}: tolerance level indicating when to stop the iteration.
\item \texttt{rcovformula}: formula to calculate the estimated residual covariance
matrix (see details). Default is equal to 1.
\item \texttt{formula3sls}: formula for calculating the 3sls estimator, one of ``GLS'',
``IV'', ``GMM'', ``Schmidt'', or ``Eviews'' (see details.)
\item \texttt{probdfsys}: use the degrees of freedom of the whole system
(in place of the degrees of freedom of the single equation to calculate probability
values for the t-test of individual parameters.
\item \texttt{single.eq.sigma}: use different $\sigma^2$ for each single
equation to calculate the covariance matrix and the standard errors of the coefficients.
\item \texttt{solvetol}: tolerance level for detecting linear dependencies when
inverting a matrix or calculating a determinant. Default is \texttt {solvetol}=.Machine\$double.eps.
\item \texttt{saveMemory}: logical. Save memory by omitting some calculation that are
not crucial for the basic estimate (e.g McElroy's $R^2$).
\end{itemize}
\subsubsection{Details}
The matrix \texttt{TX} transforms the regressor matrix ($X$) by $X\ast=X \times TX$. Thus,
the vector of coefficients is now $b=TX \times b\ast$ where $b$ is the original(stacked)
vector of all coefficients and $b\ast$ is the new coefficient vector that is estimated instead.
Thus, the elements of vector $b$ and $b_i = \sum_j TX_{ij}\times b_j\ast$. The $TX$ matrix can be
used to change the order of the coefficients and also to restrict coefficients (if $TX$ has
less columns than it has rows).
If iterated (with \texttt{maxit}>1), the covergence criterion is
\begin{eqnarray*}
\sqrt{\frac{\sum_i(b_{i,g}-b_{i,g-1})^2}{\sum_ib_{i,g-1}^2}} < tol
\end{eqnarray*}
where $b_{i,g}$ is the ith coefficient of the gth iteration step.
The formula (\texttt{rcovformula} to calculate the estimated covariance matrix of the residuals($\hat{\Sigma}$)can be one
of the following (see Judge et al., 1955, p.469):
if \texttt{rcovformula}= 0:
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T}
\end{eqnarray*}
if \texttt{rcovformula}= 1 or \texttt{rcovformula}='geomean':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{\sqrt{(T-k_i)\times (T-k_j)}}
\end{eqnarray*}
if \texttt{rcovformula}= 2 or \texttt{rcovformula}='Theil':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T-k_i-k_j+tr[X_i(X_i\prime X_i)^{-1}X_i\prime X_j(X_j\prime X_j)^{-1}X_j\prime]}
\end{eqnarray*}
if \texttt{rcovformula}= 3 or \texttt{rcovformula}='max':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T-max(k_i,k_j)}
\end{eqnarray*}
If $i = j$, formula 1, 2, and 3 are equal. All these three formulas yield unbiased estimators
for the diagonal elements of the residual covariance matrix. If $i neq j$, only formula 2
yields an unbiased estimator for the residual covariance matrix, but it is not necessarily
positive semidefinit. Thus, it is doubtful whether formula 2 is really superior to formula 1
(Theil, 1971, p.322).
The formulas to calculate the 3sls estimator lead to identical results
if the same instruments are used in all equations. If different instruments
are used in the different equations, only the GMM-3sls estimator (``GMM'')
and the 3sls estimator proposed by Schmidt (1990) (``Schmidt'') are consistent, whereas
``GMM'' is efficient relative to ``Schmidt'' (see Schmidt, 1990).
\subsubsection{Examples}
\subsubsection{Model}
\subsubsection{See Also}
For information about two stage least square regression, see
\Sref{2sls} and \texttt{help(2sls)}.
For information about seemingly unrelated regression, see
\Sref{sur} and \texttt{help(sur)}.
\subsubsection{Quantities of Interest}
\subsubsection{Output Values}
The output of each Zelig command contains useful information which you may
view. For example, if you run:
\begin{verbatim}
z.out <- zelig(formula=fml, model = "3sls", data)
\end{verbatim}
\noindent then you may examine the available information in \texttt{z.out} by
using \texttt{names(z.out)}, see the draws from the posterior distribution of
the \texttt{coefficients} by using \texttt{z.out\$coefficients}, and view a default
summary of information through \texttt{summary(z.out)}. Other elements
available through the \texttt{\$} operator are listed below:
\begin{itemize}
\item \texttt{rcovest}: residual covariance matrix used for estimation.
\item \texttt{mcelr2}: McElroys R-squared value for the system.
\item \texttt{h}: matrix of all (diagonally stacked) instrumental variables.
\item \texttt{formula3sls}: formula for calculating the 3sls estimator
\end{itemize}
\input{systemfit_output}
\begin{itemize}
\item \texttt{inst*}: instruments of the ith equation.
\item \texttt{h*}: matrix of instrumental variables of the ith equation.
\end{itemize}
\subsubsection{Contributors}
The \texttt{3sls} function is adapted from \input{contributors_systemfit}
\end{document}
\ No newline at end of file
diff --git a/R/systemfit/sur.tex b/R/systemfit/sur.tex
new file mode 100644
index 0000000..acf9d4b
--- /dev/null
+++ b/R/systemfit/sur.tex
@@ -0,0 +1 @@
+\documentclass[12pt]{book}%
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{fullpage}
\usepackage{multirow}
\usepackage{hyperref}%
\usepackage{bibentry}
\begin{document}
\subsection{\texttt{sur}: Seemingly Unrelated Regression}
\label{sur}
\texttt{sur} extends ordinary least squares analysis to estimate system of linear
equations with correlated error terms. The seemingly unrelated regression model can
be viewed as a special case of generalized least squares.
\subsubsection{Syntax}
\begin{verbatim}
> fml <- list ("mu1" = Y1 ~ X1,
"mu2" = Y2 ~ X2,
"mu3" = Y3 ~ X3)
> z.out<-zelig(formula = fml, model = "2sls", data = mydata)
> x.out <- setx(z.out)
> s.out <- sim(z.out, x = x.out)
\end{verbatim}
\subsubsection{Inputs}
\texttt{sur} regression specification has at least $M$ equations
($M \ge 2$) corresponding to the dependent variables ($Y_1, Y_2, \ldots, Y_M$).
\begin{itemize}
\item \texttt{formula}:a list whose elements are formulae corresponding to
the $M$ equations and their respective dependent and explanatory variables.
For example, when there are no constraints on the coefficients:
\begin{verbatim}
> fml <- list ("mu1" = Y1 ~ X1,
"mu2" = Y2 ~ X2,
"mu3" = Y3 ~ X3)
\end{verbatim}
\texttt{"mu1"} is the label for the first equation with Y1 as the dependent variable
and X1 as the explanatory variable. Similarly \texttt{"mu2"} and \texttt{"mu3"} are the
labels for the Y2 and Y3 equations.
\item \texttt{tag}: Users can also put constraints on the coefficients by using
the special function \texttt{tag}. \texttt{tag} takes two parameters. The first
parameter is the variable whose coefficient needs to be constrained and the second
parameter is label for the constrained coefficient. Each label uniquely identifies
the constrained coefficient. For example:
\begin{verbatim}
> fml <- list ("mu1" = Y1 ~ tag(Xc,"constrain1")+ X1,
"mu2" = Y2 ~ tag(Xc,"constrain1")+ X2,
"mu3" = Y3 ~ X3)
\end{verbatim}
\end{itemize}
\subsubsection{Additional Inputs}
\texttt{sur} takes the following additional inputs for model
specifications:
\begin{itemize}
\item \texttt{TX}: an optional matrix to transform the regressor
matrix and, hence, also the coefficient vector (see details). Default is \texttt{NULL}.
\item \texttt{maxiter}: maximum number of iterations.
\item \texttt{tol}: tolerance level indicating when to stop the iteration.
\item \texttt{rcovformula}: formula to calculate the estimated residual covariance
matrix (see details). Default is equal to 1.
\item \texttt{probdfsys}: use the degrees of freedom of the whole system
(in place of the degrees of freedom of the single equation to calculate probability
values for the t-test of individual parameters.
\item \texttt{solvetol}: tolerance level for detecting linear dependencies when
inverting a matrix or calculating a determinant. Default is \texttt {solvetol}= \begin{verbatim}.Machine\$double.eps.\end{verbatim}
\item \texttt{saveMemory}: logical. Save memory by omitting some calculation that are
not crucial for the basic estimate (e.g McElroy's $R^2$).
\end{itemize}
\subsubsection{Details}
The matrix \texttt{TX} transforms the regressor matrix ($X$) by $X\ast=X \times TX$. Thus,
the vector of coefficients is now $b=TX \times b\ast$ where $b$ is the original(stacked)
vector of all coefficients and $b\ast$ is the new coefficient vector that is estimated instead.
Thus, the elements of vector $b$ and $b_i = \sum_j TX_{ij}\times b_j\ast$. The $TX$ matrix can be
used to change the order of the coefficients and also to restrict coefficients (if $TX$ has
less columns than it has rows).
If iterated (with \texttt{maxit}>1), the covergence criterion is
\begin{eqnarray*}
\sqrt{\frac{\sum_i(b_{i,g}-b_{i,g-1})^2}{\sum_ib_{i,g-1}^2}}< tol
\end{eqnarray*}
where $b_{i,g}$ is the ith coefficient of the gth iteration step.
The formula (\texttt{rcovformula} to calculate the estimated covariance matrix of the residuals($\hat{\Sigma}$)can be one
of the following (see Judge et al., 1955, p.469):
if \texttt{rcovformula}= 0:
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T}
\end{eqnarray*}
if \texttt{rcovformula}= 1 or \texttt{rcovformula}='geomean':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{\sqrt{(T-k_i)\times (T-k_j)}}
\end{eqnarray*}
if \texttt{rcovformula}= 2 or \texttt{rcovformula}='Theil':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T-k_i-k_j+tr[X_i(X_i\prime X_i)^{-1}X_i\prime X_j(X_j\prime X_j)^{-1}X_j\prime]}
\end{eqnarray*}
if \texttt{rcovformula}= 3 or \texttt{rcovformula}='max':
\begin{eqnarray*}
\hat{\sigma_{ij}}= \frac{\hat{e_i}\prime\hat{e_j}}{T-max(k_i,k_j)}
\end{eqnarray*}
If $i = j$, formula 1, 2, and 3 are equal. All these three formulas yield unbiased estimators
for the diagonal elements of the residual covariance matrix. If $i neq j$, only formula 2
yields an unbiased estimator for the residual covariance matrix, but it is not necessarily
positive semidefinit. Thus, it is doubtful whether formula 2 is really superior to formula 1
(Theil, 1971, p.322).
\subsubsection{Examples}
\subsubsection{Model}
The basic seemingly unrelated regression model assumes that for each individual
observation $i$ there are $M$ dependent variables ($Y_{ij}, j=1,\ldots,M$)
each with its own regression equation:
\begin{eqnarray*}
Y_{ij} = X_{ij}'\beta_j + \epsilon_{ij}, \quad {\rm for} \quad i=1,\ldots,N \quad {\rm and} \quad j=1,\ldots,M
\end{eqnarray*}
when $X_{ij}$ is a k-vector of explanatory variables, $\beta_j$
is the coefficients of the explanatory variables,
\begin{itemize}
\item The \emph{stochastic component} is:
\begin{eqnarray*}
\epsilon_{ij} & \sim & {\cal N}(0, \sigma_{ij})
\end{eqnarray*}
where within each $j$ equation, $epsilon_{ij}$ is identically
and independently distributed for $i=1,\ldots,M$,
\begin{eqnarray*}
{\rm Var}(\epsilon_{ij})=\sigma_j \quad {\rm and} \quad {\rm Cov}(\epsilon_{ij}, \epsilon_{i\prime j})= 0, \quad {\rm for} \quad i \neq i\prime,\quad {\rm and} \quad j=1,\ldots,M
\end{eqnarray*}
However, the error terms for the \emph{ith} observation can be correlated across equations
\begin{eqnarray*}
{\rm Cov}(\epsilon_{ij}, \epsilon_{ij\prime})\neq 0, \quad {\rm for} \quad j \neg j\prime, \quad {\rm and} \quad i=1,\ldots,N
\end{eqnarray*}
\item The \emph{systematic component} is:
\begin{eqnarray*}
\mu_{ij}= E(Y_ij)= X_{ij}\beta_j, \quad {\rm for} \quad i=1,\ldots,N, \quad {\rm and} \quad j=1,\ldots,M
\end{eqnarray*}
\end{itemize}
\subsubsection{See Also}
For information about two stage least squares regression, see
\Sref{2sls} and \texttt{help(2sls)}.
For information about three stage least squares regression, see
\Sref{3sls} and \texttt{help(3sls)}.
\subsubsection{Quantities of Interest}
\subsubsection{Output Values}
The output of each Zelig command contains useful information which you may
view. For example, if you run:
\begin{verbatim}
z.out <- zelig(formula=fml, model = "sur", data)
\end{verbatim}
\noindent then you may examine the available information in \texttt{z.out} by
using \texttt{names(z.out)}, see the draws from the posterior distribution of
the \texttt{coefficients} by using \texttt{z.out\$coefficients}, and view a default
summary of information through \texttt{summary(z.out)}. Other elements
available through the \texttt{\$} operator are listed below:
\begin{itemize}
\item \texttt{rcovest}: residual covariance matrix used for estimation.
\item \texttt{mcelr2}: McElroys R-squared value for the system.
\end{itemize}
\input{systemfit_output}
\begin{itemize}
\item \texttt{maxiter}: maximum number of iterations.
\item \texttt{tol}: tolerance level indicating when to stop the iteration.
\end{itemize}
\subsubsection{Contributors}
The \texttt{sur} function is adapted from \input{contributors_systemfit}
\end{document}
\ No newline at end of file
diff --git a/R/terms.multiple.R b/R/terms.multiple.R
new file mode 100644
index 0000000..c52b6a4
--- /dev/null
+++ b/R/terms.multiple.R
@@ -0,0 +1,98 @@
+terms.multiple<-function(x, data=NULL,...){
+ object <- x
+ termsexist<-attr(object,"terms")
+ if(!(is.null(termsexist)))
+ return (termsexist)
+
+ nreq<-nrConstr<-nrEquationsNew<-0
+ constr<-XconstrEqn<-variables<-termlabels<-depVars<-objectNew<-intercAttr<-depFactors<-list()
+ depFactorVar<-depLevels<-namesConstr<-c()
+ namesOfEquations<- names(object)
+ nrEquations <-length(object)
+ "%w/o%" <- function(x,y) x[!x %in% y]
+
+ for (i in 1:nrEquations){
+ TT<-terms.formula(object[[i]], specials=c("id","tag"))
+ eqni<-object[[i]]
+ namei<-namesOfEquations[[i]]
+ tagattr<-attr(TT,"specials")$tag
+ hastag<-!(is.null(tagattr))
+ if (hastag){
+ constrTmp<-c()
+ for(j in 1:length(tagattr)){
+ nrConstr<-nrConstr+1
+ if(length(eqni)==3)
+ ind<-tagattr[[j]]-1
+ else
+ ind<-tagattr[[j]]
+ vind<-tagattr[[j]]+1
+ parsedTag<-attr(TT,"variables")[[vind]]
+ varname<-as.character(parsedTag[[2]])
+ label<-as.character(parsedTag[[3]])
+ namesConstr<-c(namesConstr,label)
+ constr[[nrConstr]]<-c(i,label,varname)
+ attr(TT,"term.labels")[[ind]]<-varname
+ attr(TT,"variables")[[vind]]<-as.name(varname)
+ constrTmp<-c(constrTmp,varname)
+ }
+ XconstrEqn[[i]]<-constrTmp
+ }
+ if (length(eqni)==3){
+ nrEquationsNew<-nrEquationsNew+1
+ objectNew[[namei]]<-eqni
+ nreq=nreq+1
+ lhs<-eqni[[2]]
+ if (length(lhs)>1 && lhs[[1]]=="id"){
+ depVars[[namei]]<-lhs[[3]]
+ depFactorVar<-c(depFactors,deparse(lhs[[2]]))
+ depLevels<-c(depLevels,lhs[[3]])
+ }else
+ depVars[[namei]]<-deparse(eqni[[2]])
+ }else{
+ nrEquationsNew<-nrEquationsNew+1
+ objectNew[[namei]]<-eqni
+ }
+ variables[[namei]]<-attr(TT,"variables")
+ termlabels[[namei]]<-attr(TT,"term.labels")
+ intercAttr[[namei]]<-attr(TT,"intercept")
+ }
+ namesOfEquations<-names(objectNew)
+ myattr<-list()
+ result<-objectNew
+ if(length(constr)>0){
+ namesConstr<-unique(namesConstr)
+ constraints<-matrix(NA,nrow=nrEquationsNew,ncol=length(namesConstr),dimnames=list(namesOfEquations,namesConstr))
+ for(i in 1:length(constr)){
+ constri<-constr[[i]]
+ eqind<-constri[[1]]
+ eq<-namesOfEquations[as.numeric(eqind)]
+ lab<-constri[[2]]
+ constraints[eq,lab]<-constri[[3]]
+ }
+ }else
+ constraints<-FALSE
+
+ indVars<-unique(unlist(termlabels))
+ if(length(depFactorVar) !=0)
+ depFactors<-list("depFactorVar"=unique(unlist(depFactorVar)),"depLevels"=depLevels)
+ else
+ depFactors<-FALSE
+
+ whiche<-which(lapply(termlabels,length)!=0)
+ myattr$systEqns<-names(whiche)
+ myattr$ancilEqns<-"%w/o%"(namesOfEquations,myattr$systEqns)
+
+ myattr$variables<-variables
+ myattr$term.labels<-termlabels
+ myattr$indVars<-indVars
+
+ myattr$depVars<-depVars
+ myattr$depFactors<-depFactors
+ myattr$constraints<-constraints
+ myattr$response<-1
+ myattr$intercept<-intercAttr
+ attributes(result)<-myattr
+ names(result)<-namesOfEquations
+ class(result)<-c("terms","multiple","list")
+ return(result)
+}
diff --git a/R/ternaryplot.R b/R/ternaryplot.R
new file mode 100644
index 0000000..7ffc497
--- /dev/null
+++ b/R/ternaryplot.R
@@ -0,0 +1,89 @@
+ternaryplot <- function (x, scale = 1, dimnames = NULL, dimnames.position = c("corner",
+ "edge", "none"), dimnames.color = "black", id = NULL, id.color = "black",
+ coordinates = FALSE, grid = TRUE, grid.color = "gray", labels = c("inside",
+ "outside", "none"), labels.color = "darkgray", border = "black",
+ bg = "white", pch = 19, cex = 1, prop.size = FALSE, col = "red",
+ main = "ternary plot", ...) {
+
+ ## From vcd (Version 0.1-3.3). Function by David Meyer
+ labels <- match.arg(labels)
+ if (grid == TRUE)
+ grid <- "dotted"
+ if (coordinates)
+ id <- paste("(", round(x[, 1] * scale, 1), ",", round(x[,
+ 2] * scale, 1), ",", round(x[, 3] * scale, 1), ")",
+ sep = "")
+ dimnames.position <- match.arg(dimnames.position)
+ if (is.null(dimnames) && dimnames.position != "none")
+ dimnames <- colnames(x)
+ if (is.logical(prop.size) && prop.size)
+ prop.size <- 3
+ if (ncol(x) != 3)
+ stop("Need a matrix with 3 columns")
+ if (any(x) < 0)
+ stop("X must be non-negative")
+ s <- rowSums(x)
+ if (any(s <= 0))
+ stop("each row of X must have a positive sum")
+ x <- x/s
+ top <- sqrt(3)/2
+ par(plt = c(0.06, 0.94, 0.15, 0.87))
+ plot.new()
+ xlim <- c(-0.03, 1.03)
+ ylim <- c(0, top)
+ par(usr = c(xlim, ylim), oma = c(0, 0, 1, 0))
+ plot.window(xlim = xlim, ylim = ylim, asp = 1)
+ eps <- 0.01
+ polygon(c(0, 0.5, 1), c(0, top, 0), col = bg, xpd = NA, border = border,
+ ...)
+ title(main, outer = TRUE, line = -1)
+ if (dimnames.position == "corner") {
+ axis(1, at = c(-0.03, 1.03), labels = dimnames[1:2],
+ tick = FALSE, font = 2)
+ axis(3, at = 0.5, labels = dimnames[3], tick = FALSE,
+ font = 2)
+ }
+ if (dimnames.position == "edge") {
+ shift <- eps * if (labels == "outside")
+ 8
+ else 0
+ text(0.25 - 2 * eps - shift, 0.5 * top + shift, dimnames[2],
+ srt = 60, col = dimnames.color)
+ text(0.75 + 3 * eps + shift, 0.5 * top + shift, dimnames[1],
+ srt = -60, col = dimnames.color)
+ text(0.5, 0, dimnames[3], pos = 1, offset = 0.5 + 30 *
+ shift, xpd = NA, col = dimnames.color)
+ }
+ if (is.character(grid))
+ for (i in 1:4 * 0.2) {
+ lines(c(1 - i, (1 - i)/2), c(0, 1 - i) * top, lty = grid,
+ col = grid.color)
+ lines(c(1 - i, 1 - i + i/2), c(0, i) * top, lty = grid,
+ col = grid.color)
+ lines(c(i/2, 1 - i + i/2), c(i, i) * top, lty = grid,
+ col = grid.color)
+ if (labels == "inside") {
+ text((1 - i) * 3/4 - eps, (1 - i)/2 * top, i *
+ scale, col = labels.color, srt = 120)
+ text(1 - i + i/4 + eps, i/2 * top - eps, (1 -
+ i) * scale, col = labels.color, srt = -120)
+ text(0.5, i * top + eps, i * scale, col = labels.color)
+ }
+ if (labels == "outside") {
+ text((1 - i)/2 - 6 * eps, (1 - i) * top, (1 -
+ i) * scale, col = labels.color)
+ text(1 - (1 - i)/2 + 3 * eps, (1 - i) * top +
+ 5 * eps, i * scale, srt = -120, col = labels.color)
+ text(i + eps, 0, (1 - i) * scale, pos = 1, offset = 1.5,
+ srt = 120, xpd = NA, col = labels.color)
+ }
+ }
+ xp <- x[, 2] + x[, 3]/2
+ yp <- x[, 3] * top
+ points(xp, yp, pch = pch, col = col, cex = if (prop.size)
+ prop.size * (s/max(s))
+ else cex, ...)
+ if (!is.null(id))
+ text(xp, yp, as.character(id), pos = 1, offset = 0.5 *
+ cex, col = id.color)
+}
diff --git a/R/vdc.R b/R/vdc.R
new file mode 100644
index 0000000..58d370f
--- /dev/null
+++ b/R/vdc.R
@@ -0,0 +1,244 @@
+zeligListModels<-function(inZeligOnly=T) {
+ if (inZeligOnly) {
+ tmp = ls(envir=asNamespace("Zelig"),pattern="^zelig2")
+ } else {
+ tmp = c( ls(envir=asNamespace("Zelig"),pattern="^zelig2"),
+ apropos("zelig2",mode="function"))
+ }
+ sub("zelig2","", tmp)
+}
+
+zeligInstalledModels<-function(inZeligOnly=T,schemaVersion="1.1") {
+ chkpkgs<-function(name) {
+ zd=zeligDescribeModelXML(name,schemaVersion=schemaVersion)
+ if (is.null(zd)) {
+ return (FALSE)
+ }
+ zdpd= zeligModelDependency(name)[,1]
+ if (is.null(zdpd)) {
+ return(TRUE)
+ }
+ ow=options(warn=-1)
+ ret = (class(try(sapply(zdpd,function(x)require(x,character.only=T)),silent=T))
+ !="try-error")
+ options(ow)
+ return (ret)
+ }
+ models<-zeligListModels(inZeligOnly=inZeligOnly)
+ # Not being the trusting sort, lets check to see if we can run
+ # a dummy formula. If not -- almost always means that something
+ # required() is missing
+ tmpModels<-sapply(models,chkpkgs)
+ models[which(tmpModels)]
+}
+
+zeligGetSpecial<-function(modelName) {
+ modelDesc = zeligDescribeModel(modelName)
+ return(modelDesc$parameters[[1]]$specialFunction)
+}
+
+zeligModelDependency<-function(modelName,repos="") {
+ zd= zeligDescribeModel(modelName)
+
+ if (is.null(zd)) { return (NULL) }
+
+ zdpd=zd[which(names(zd)=="package")]
+
+ cbind(sapply(zdpd,function(x)x$name),
+ sapply(zdpd,function(x){if (is.null(x$CRAN))
+ {rv<-repos} else{rv<-x$CRAN};rv;}))
+ }
+
+
+zeligDescribeModel<-function(name,force=F,schemaVersion="1.1") {
+ res=try(eval(call(paste("describe.",name,sep=""))),silent=T)
+ if (inherits(res,"try-error")) {
+ if (force) {
+ res=describe.default()
+ } else {
+ res=NULL
+ }
+ }
+# res$name<-name # only here we have access to the model name, so add it to the list.
+ if(!is.null(res)){
+ res<-check.full(res,name)
+ }
+ return(res)
+}
+
+zeligDescribeModelXML<-function(modelName,force=F,schemaVersion="1.1") {
+ zd = zeligDescribeModel(modelName,force,schemaVersion)
+ if (is.null(zd)) {
+ return(NULL)
+ } else {
+ return(zmodel2string(zd))
+ }
+
+}
+
+
+zmodel2string<-function(x) {
+ xmlList(x)
+}
+printZeligSchemaInstance<-function(filename=NULL, serverName=NULL,vdcAbsDirPrefix=NULL){
+ # open connection
+ schemaURL<-'http://gking.harvard.edu/zelig';
+ if (is.null(serverName)) {
+ serverName<-system('hostname -f', intern=T)
+ }
+ if (is.null(vdcAbsDirPrefix)){
+ locationURL<-paste('http://', serverName, '/VDC/schema/analysis/ZeligInterfaceDefinition.xsd',sep="");
+ } else {
+ locationURL<-paste('file://', vdcAbsDirPrefix, '/VDC/schema/analysis/ZeligInterfaceDefinition.xsd',sep="");
+ }
+ schemaLocation<-paste(schemaURL, ' ', locationURL, sep='');
+ con<-"";
+ if (!is.null(filename)){
+ con<-file(filename,"w");
+ }
+ cat(file=con, "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n<zelig xmlns=\"",schemaURL,"\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" xsi:schemaLocation=\"",schemaLocation,"\">", sep="");
+ mssg<- sapply(zeligInstalledModels(),function(x){cat(file=con,zmodel2string(zeligDescribeModel(x)),sep="")},simplify=F);
+ cat(file=con,"\n</zelig>\n",sep="");
+}
+
+
+xmlList<-function(z){
+ if(is.null(z))
+ return ("")
+ mins<-c()
+ maxs<-c()
+ for(i in 1:length(z$parameters)){
+ mins<-c(mins,z$parameters[[i]]$equations[[1]])
+ maxs<-c(maxs,z$parameters[[i]]$equations[[2]])
+ }
+ min<-sum(mins)
+ if(any(!is.finite(maxs)))
+ max<-Inf
+ else
+ max<-sum(maxs)
+ if(max !=1)
+ return("")
+ res<-paste("<model name=",'"',z$name,'"'," label=",'"',categories()[[z$category]],'"', sep="")
+ if(!(is.na(z$parameters[[1]]$specialFunction)))
+ res<-paste(res," specialFunction=",'"',z$parameters[[1]]$specialFunction,'"',sep="")
+ res<-paste(res,">\n",sep="")
+ res<-paste(res,"<description>",z$description, "</description>\n",sep="")
+ if(z$name=="irtkd")
+ res<-paste(res,"<helpLink url=","http://gking.harvard.edu/zelig/docs/_TT_irtkd_TT__tex2htm.html",sep="")
+ else
+ res<-paste(res,"<helpLink url=",'"',modelURL(z$name,z$description),'"',sep="")
+
+ res<-paste(res,"/>\n", sep="")
+ if(any(!(is.null(z$packageDependency)))){
+ res<-paste(res,"<packageDependency",sep="")
+ if(!(is.na(z$packageDependency$name)))
+ res<-paste(res," name= ",'"',z$packageDependency$name,'"',sep="")
+ if(!(is.na(z$packageDependency$version)))
+ res<-paste(res," version= ",'"',z$packageDependency$version,'"',sep="")
+ if(!(is.na(z$packageDependency$relationship)))
+ res<-paste(res," relationship= ",'"',z$packageDependency$relationship,'"',sep="")
+ if(!(is.na(z$packageDependency$CRAN)))
+ res<-paste(res," CRAN= ",'"',z$packageDependency$CRAN,'"',sep="")
+ res<-paste(res,"/>\n",sep="")
+ }
+ res<-paste(res,"<formula minEquations=",'"',min,'"',sep="")
+ if(is.finite(max))
+ res<-paste(res," maxEquations=",'"',max,'"',sep="")
+ if(max==1)
+ res<-paste(res," simulEq=",'"',0,'"',sep="")
+ res<-paste(res,"/>\n",sep="")
+
+ res<-paste(res,"<equation name=",'"',names(z$parameters)[[1]],'"',"/>\n",sep="")
+ if(z$parameters[[1]]$depVar){
+ res<-paste(res,"<outcome",sep="")
+ if(!is.na(z$parameters[[1]]$specialFunction))
+ {
+ if(is.finite(z$parameters[[1]]$varInSpecialFunction[[2]] ))
+ res<-paste(res," maxVar=",'"',z$parameters[[1]]$varInSpecialFunction[[2]],'"',sep="")
+ res<-paste(res," minVar=",'"',z$parameters[[1]]$varInSpecialFunction[[1]],'"',sep="")
+ }
+ else
+ {
+ if(z$parameters[[1]]$depVar == TRUE)
+ res<-paste(res," minVar=",'"',1,'"',sep="")
+ else
+ res<-paste(res," minVar=",'"',0,'"'," maxVar=",'"',0,'"',sep="")
+ }
+ res<-paste(res,">\n")
+
+
+ for(i in 1:length(modeling.types()$depVar[[z$category]]))
+ res<-paste(res,"<modelingType>",modeling.types()$depVar[[z$category]][[i]],"</modelingType>\n",sep="")
+
+ res<-paste(res,"</outcome>\n",sep="")
+ }
+ #explanatory
+ if(z$parameters[[1]]$expVar){
+ res<-paste(res,"<explanatory ")
+ if(z$parameters[[1]]$expVar == TRUE)
+ res<-paste(res," minVar=",'"',1,'"',sep="")
+ else
+ res<-paste(res," minVar=",'"',0,'"'," maxVar=",'"',0,'"',sep="")
+ res<-paste(res,">\n")
+
+
+ for(i in 1:length(modeling.types()$expVar[[z$category]]))
+ res<-paste(res,"<modelingType>", modeling.types()$expVar[[z$category]][[i]],"</modelingType>\n",sep="")
+ res<-paste(res,"</explanatory>\n",sep="")
+ }
+ res<-paste(res,"</equation>\n",sep="")
+ res<-paste(res,"</formula>\n",sep="")
+ res<-paste(res,"<setx maxSetx=",'"',2,'"',"/>\n",sep="")
+ res<-paste(res,"</model>\n")
+ return(res)
+}
+
+
+
+check.full<-function(z,name){
+ # we suppose that describe.model pass the check
+ z$name<-name
+ if(is.null(z$package))
+ z$package<-list(name=NA,version=NA, CRAN=NA)
+
+ for (i in length(z$parameters)){
+ if(is.null(z$parameters[[i]]$specialFunction)) z$parameters[[i]]$specialFunction<-NA
+ if(is.null(z$parameters[[i]]$varInSpecial)) z$parameters[[i]]$varInSpecial<-NA
+}
+ return(z)
+
+}
+
+
+
+
+
+modelURL<-function(modelName,modelDesc){
+ baseUrl<-"http://gking.harvard.edu/zelig/docs/"
+ spec<-"_TT_"
+ res<-paste(baseUrl,spec,modelName,spec,"_",substr(modelDesc,0, 13-nchar(modelName)) ,".html",sep="")
+ res<-gsub(".","_",res,fixed=TRUE)
+ res<-gsub(" ","_",res,fixed=TRUE)
+ res
+}
+
+modeling.types <-function(){
+ res<-list(
+ expVar=list(continuous=c("continuous","discrete","nominal","ordinal","binary"),
+ dichotomous=c("continuous","discrete","nominal","ordinal","binary"),
+ ordinal=c("continuous","discrete","nominal","ordinal","binary"),
+ bivariate.dichotomous=c("continuous","discrete","nominal","ordinal","binary"),
+ multinomial=c("continuous","discrete","nominal","ordinal","binary"),
+ event.count=c("continuous","discrete","nominal","ordinal","binary"),
+ censored=c("continuous","discrete","nominal","ordinal","binary")),
+ depVar=list(
+ continuous="continous",
+ dichotomous="binary",
+ ordinal="ordinal",
+ bivariate.dichotomous="binary",
+ multinomial=c("nominal","ordinal"),
+ event.count="discrete",
+ censored="continuous"))
+ res
+}
+
diff --git a/R/zelig.R b/R/zelig.R
index e1140b1..e88b142 100644
--- a/R/zelig.R
+++ b/R/zelig.R
@@ -1,72 +1,2 @@
-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 <- zelig.call <- match.call(expand.dots = TRUE)
- if (missing(by))
- by <- NULL
- N <- M <- 1
- object <- list()
- if (!is.data.frame(data))
- M <- length(data)
- if (M > 1)
- dat <- data[[1]]
- else
- dat <- data
- if (!is.null(by)) {
- if (any(as.character(by) %in% c(formula[[2]], formula[[3]])))
- stop("the variable selected for subsetting cannot be called in the formula.")
- idx <- dat[,by]
- mf$by <- NULL
- lev <- sort(unique(idx))
- N <- length(lev)
- }
- mf <- do.call(fn1, list(formula, model, dat, N, ...))
- for (i in 1:N) {
- if (N > 1) {
- dat <- list()
- if (M > 1) {
- for (j in 1:M)
- dat[[j]] <- data[[j]][idx == lev[i],]
- }
- else
- dat <- data[idx == lev[i],]
- }
- else
- dat <- data
- obj <- list()
- for (j in 1:M) {
- if (M > 1)
- d <- dat[[j]]
- else
- d <- dat
- if (is.data.frame(d)) {
- d <- d[complete.cases(model.frame(as.formula(formula), d)),]
- mf$data <- d
- res <- eval(as.call(mf))
- if (exists(fn2))
- res <- do.call(fn2, list(res = res, fcall = mf,
- zcall = as.list(zelig.call)))
- res$call <- zelig.call
- res$data <- res$call$data
- res$zelig <- model
- if (M > 1)
- obj[[j]] <- res
- else
- obj <- res
- }
- }
- if (M > 1)
- class(obj) <- "MI"
- if (N > 1)
- object[[i]] <- obj
- else
- object <- obj
- }
- if (N > 1) {
- class(object) <- "strata"
- names(object) <- lev
- }
- object
-}
+zelig <- function(formula, ...)
+ UseMethod("zelig")
diff --git a/R/zelig.R b/R/zelig.default.R
similarity index 81%
copy from R/zelig.R
copy to R/zelig.default.R
index e1140b1..5f59898 100644
--- a/R/zelig.R
+++ b/R/zelig.default.R
@@ -1,9 +1,11 @@
-zelig <- function(formula, model, data, by = NULL, ...) {
+zelig.default <- function(formula, model, data, by = NULL, save.data =
+ FALSE, ...) {
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 <- zelig.call <- match.call(expand.dots = TRUE)
+ zelig.call[[1]] <- as.name("zelig")
if (missing(by))
by <- NULL
N <- M <- 1
@@ -42,14 +44,18 @@ zelig <- function(formula, model, data, by = NULL, ...) {
else
d <- dat
if (is.data.frame(d)) {
- d <- d[complete.cases(model.frame(as.formula(formula), d)),]
+ d <- d[complete.cases(model.frame(mf$formula, data=d,
+ na.action = na.pass)),]
mf$data <- d
res <- eval(as.call(mf))
if (exists(fn2))
res <- do.call(fn2, list(res = res, fcall = mf,
zcall = as.list(zelig.call)))
res$call <- zelig.call
- res$data <- res$call$data
+ if (save.data)
+ res$data <- d
+ else
+ res$data <- res$call$data
res$zelig <- model
if (M > 1)
obj[[j]] <- res
diff --git a/R/zelig2blogit.R b/R/zelig2blogit.R
index 63fe184..c9ed318 100644
--- a/R/zelig2blogit.R
+++ b/R/zelig2blogit.R
@@ -1,5 +1,4 @@
-zelig2blogit <- function(formula, model, data, M, constrain = NULL,
- omit = NULL, constant = 3, ...) {
+zelig2blogit <- function(formula, model, data, M, constant = 3, ...) {
check <- library()
if(any(check$results[,"Package"] == "VGAM"))
require(VGAM)
@@ -9,11 +8,12 @@ zelig2blogit <- function(formula, model, data, M, constrain = NULL,
mf[[1]] <- VGAM::vglm
mf$family <- as.name("blogit")
mf$... <- NULL
- tmp <- cmvglm(formula, model, constrain, omit, constant, 3)
+ formula<-parse.formula(formula,model)
+ tmp <- cmvglm(formula, model, constant, 3)
mf$formula <- tmp$formula
mf$constraints <- tmp$constraints
blogit <<- function() binom2.or(zero=NULL)
- mf$model <- mf$constrain <- mf$omit <- mf$constant <- mf$M <- NULL
+ mf$model <- mf$constant <- mf$M <- NULL
mf$robust <- NULL
as.call(mf)
}
diff --git a/R/zelig2bprobit.R b/R/zelig2bprobit.R
index 2c28116..10bcaa3 100644
--- a/R/zelig2bprobit.R
+++ b/R/zelig2bprobit.R
@@ -1,5 +1,4 @@
-zelig2bprobit <- function(formula, model, data, M, constrain = NULL,
- omit = NULL, constant = 3, ...) {
+zelig2bprobit <- function(formula, model, data, M, constant = 3, ...) {
check <- library()
if(any(check$results[,"Package"] == "VGAM"))
require(VGAM)
@@ -8,10 +7,12 @@ zelig2bprobit <- function(formula, model, data, M, constrain = NULL,
mf <- match.call(expand.dots = TRUE)
mf[[1]] <- VGAM::vglm
mf$family <- as.name("bprobit")
- tmp <- cmvglm(formula, model, constrain, omit, constant, 3)
- mf$formula <- tmp$formula
+ formula<-parse.formula(formula,model)
+ tmp <- cmvglm(formula, model, constant, 3)
+ mf$formula <- tmp$formula
mf$constraints <- tmp$constraints
bprobit <<- function() binom2.rho(zero=NULL)
- mf$model <- mf$constrain <- mf$omit <- mf$constant <- mf$M <- NULL
+ mf$model <- mf$constant <- mf$M <- NULL
+
as.call(mf)
}
diff --git a/R/zelig2gamma.R b/R/zelig2gamma.R
index afd3ff6..ab08b2d 100644
--- a/R/zelig2gamma.R
+++ b/R/zelig2gamma.R
@@ -1,6 +1,7 @@
zelig2gamma <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- mf$robust <- NULL
+ mf$M <- mf$robust <- NULL
+ mf$model <- FALSE
mf[[1]] <- glm
mf$family <- Gamma
as.call(mf)
diff --git a/R/zelig2logit.R b/R/zelig2logit.R
index 8585cb8..9c32781 100644
--- a/R/zelig2logit.R
+++ b/R/zelig2logit.R
@@ -1,8 +1,7 @@
zelig2logit <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- mf$robust <- NULL
- mf$formula <- as.formula(paste("cbind(", formula[[2]], ", 1 -", formula[[2]], ")",
- "~", deparse(formula[[3]], width.cutoff=500), sep = ""))
+ mf$M <- mf$robust <- NULL
+ mf$model <- FALSE
mf[[1]] <- stats::glm
mf$family <- binomial(link="logit")
as.call(mf)
diff --git a/R/zelig2ls.R b/R/zelig2ls.R
index e2af63f..c18179f 100644
--- a/R/zelig2ls.R
+++ b/R/zelig2ls.R
@@ -1,6 +1,7 @@
zelig2ls <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- mf$robust <- NULL
+ mf$model <- FALSE
+ mf$M <- mf$robust <- NULL
mf[[1]] <- lm
as.call(mf)
}
diff --git a/R/zelig2mlogit.R b/R/zelig2mlogit.R
index 871e811..c603a3b 100644
--- a/R/zelig2mlogit.R
+++ b/R/zelig2mlogit.R
@@ -7,17 +7,13 @@ zelig2mlogit <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
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) {
- ntmp <- array()
- for (i in 1:mf$M)
- ntmp[i] <- length(unique(na.omit((eval(mf$formula[[2]], mf$data[[i]]))))) $
- ndim <- max(ntmp)
- }
- tmp <- cmvglm(mf$formula, mf$model, mf$equal, mf$zeros, mf$ones, ndim)
+ formula<-parse.formula(formula,model,data)
+ tt<-terms(formula)
+ fact<-attr(tt,"depFactors")$depFactorVar
+ ndim<-length(attr(tt,"depFactors")$depLevels)
+ tmp <- cmvglm(formula, mf$model, mf$ones, ndim,data,fact)
mf$formula <- tmp$formula
mf$constraints <- tmp$constraints
- mf$model <- mf$equal <- mf$ones <- mf$zeros <- mf$M <- NULL
+ mf$model <- mf$M <- NULL
as.call(mf)
}
diff --git a/R/zelig2normal.R b/R/zelig2normal.R
index 5a9f884..8c97a06 100644
--- a/R/zelig2normal.R
+++ b/R/zelig2normal.R
@@ -1,6 +1,7 @@
zelig2normal <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- mf$robust <- NULL
+ mf$M <- mf$robust <- NULL
+ mf$model <- FALSE
mf[[1]] <- glm
mf$family <- gaussian
as.call(mf)
diff --git a/R/zelig2poisson.R b/R/zelig2poisson.R
index 11b4e2e..9d967ac 100644
--- a/R/zelig2poisson.R
+++ b/R/zelig2poisson.R
@@ -1,6 +1,7 @@
zelig2poisson <- function(formula, model, data, M, ...) {
mf <- match.call(expand.dots = TRUE)
- mf$model <- mf$M <- mf$robust <- NULL
+ mf$M <- mf$robust <- NULL
+ mf$model <- FALSE
mf[[1]] <- glm
mf$family <- poisson
as.call(mf)
diff --git a/R/zelig2probit.R b/R/zelig2probit.R
index e9ddc53..f3f40d3 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 <- mf$probit <- NULL
- mf$formula <- as.formula(paste("cbind(", formula[[2]], ", 1 -", formula[[2]], ")",
- "~", deparse(formula[[3]], width.cutoff=500), sep = ""))
+ mf$M <- mf$probit <- NULL
+ mf$model <- FALSE
mf[[1]] <- stats::glm
mf$family <- binomial(link="probit")
as.call(mf)
diff --git a/R/zelig2tobit.R b/R/zelig2tobit.R
new file mode 100644
index 0000000..80dbf27
--- /dev/null
+++ b/R/zelig2tobit.R
@@ -0,0 +1,59 @@
+zelig2tobit <- function(formula, model, data, M, ...) {
+ mf <- match.call(expand.dots = TRUE)
+ mf$model <- mf$M <- NULL
+ require(survival)
+ mf[[1]] <- survival::survreg
+ mf$dist <- "gaussian"
+ if (is.null(mf$above))
+ above <- Inf
+ else {
+ above <- mf$above
+ mf$above <- NULL
+ }
+ if (is.null(mf$below))
+ below <- 0
+ else {
+ below <- mf$below
+ mf$below <- NULL
+ }
+
+ ## Fixing the call for robust SEs
+ 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]], width.cutoff=500)),
+ 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]], width.cutoff=500)),
+ paste("+", " cluster(1:nrow(",
+ deparse(formula[[2]]),"))")))
+
+ ## Fixing the call for left censoring
+ if (length(grep("Surv", as.character(formula[[2]]))) == 0) {
+ if (above == Inf & is.numeric(below) & below != -Inf) {
+ tt <- "left"
+ }
+ else if (below == -Inf & above == Inf) {
+ stop("No censoring, use zelig(..., model = \"normal\") instead!")
+ }
+ else if (below == -Inf & is.numeric(above) & above != Inf) {
+ stop("Right censored data not supported in a tobit model.")
+ }
+ else if (is.numeric(above) & is.numeric(below) & above != Inf & below != -Inf) {
+ stop("Interval censored data not suppored in a tobit model.")
+ }
+ mf$formula[[2]] <- call("Surv", formula[[2]],
+ call("<", below, formula[[2]]),
+ type = "left")
+ }
+ as.call(mf)
+}
diff --git a/R/zelig3glm.R b/R/zelig3glm.R
index 3349769..a29b4a2 100644
--- a/R/zelig3glm.R
+++ b/R/zelig3glm.R
@@ -1,7 +1,7 @@
zelig3logit <- zelig3probit <- zelig3normal <- zelig3gamma <-
zelig3poisson <- zelig3negbin <- zelig3relogit <-
function(res, fcall = NULL, zcall = NULL) {
- rob <- zcall$robust
+ rob <- eval(zcall$robust)
if (!is.null(rob)) {
require(sandwich)
if (is.list(rob)) {
diff --git a/R/zelig3relogit.R b/R/zelig3relogit.R
index 3d9a482..2528fa5 100644
--- a/R/zelig3relogit.R
+++ b/R/zelig3relogit.R
@@ -1,27 +1,34 @@
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)
+ obj <- list()
+ obj$lower.estimate <- zelig3relogit(res$lower.estimate, fcall =
+ fcall, zcall = zcall)
+ obj$upper.estimate <- zelig3relogit(res$upper.estimate, fcall =
+ fcall, zcall = zcall)
+ obj$upper.estimate$call <- obj$lower.estimate$call <- as.call(zcall)
+ class(obj) <- class(res)
+ return(obj)
}
+
+ zcall$robust <- eval(zcall$robust)
if (is.null(zcall$robust)) {
- if (is.null(zcall$weighting))
- rob <- FALSE
- else if (zcall$weighting) {
+ if (res$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) {
+ else if (is.logical(zcall$robust)) {
+ if (!zcall$robust & res$weighting) {
rob <- TRUE
warning("robust is set to TRUE because weighting is used")
}
+ else
+ rob <- zcall$robust
+ }
else
rob <- zcall$robust
if (is.list(rob)) {
@@ -29,24 +36,15 @@ zelig3relogit <- function(res, fcall = NULL, zcall = NULL) {
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
- }
+ 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")
+ class(res) <- c("relogit", "glm.robust")
}
return(res)
}
diff --git a/README b/README
index 3341c47..b562ee3 100644
--- a/README
+++ b/README
@@ -1,3 +1,24 @@
+2.5-2 (February 3, 2006): Stable release for R 2.0.0-2.2.1. Fixed bugs related to
+ VDC GUI. First level dependencies are the same as in version 2.5-1.
+
+2.5-1 (January 31, 2006): Stable release for R 2.0.0-2.2.1. Added methods for
+ multiple equation models. Fixed bugs related to robust estimation and upgrade of
+ sandwich and zoo packages. Revised setx() to use environments. Added
+ current.packages() to retrieve version of packages upon which Zelig depends.
+ First level version dependencies are as follows:
+ MASS 7.2-24
+ boot 1.2-24
+ VGAM 0.6-7
+ mvtnorm 0.7-2
+ survival 2.20
+ sandwich 1.1-0
+ zoo 1.0-4
+ MCMCpack 0.6-6
+ coda 0.10-3
+
+2.4-6 (October 27, 2005): Stable release for R 2.0.0-2.2.0. Fixed bug
+ related to simulation for Bayesian Normal regression.
+
2.4-5 (October 18, 2005): Stable release for R 2.0.0-2.2.0. Updated
installation instructions.
diff --git a/data/Zelig.url.tab b/data/Zelig.url.tab
index 2a72c0b..ac405dc 100644
--- a/data/Zelig.url.tab
+++ b/data/Zelig.url.tab
@@ -43,4 +43,6 @@ poisson.bayes http://gking.harvard.edu/zelig/docs/_TT_poisson_bayes_TT.html
probit http://gking.harvard.edu/zelig/docs/_TT_probit_TT__Probit.html
probit.bayes http://gking.harvard.edu/zelig/docs/_TT_probit_bayes_TT__B.html
relogit http://gking.harvard.edu/zelig/docs/_TT_relogit_TT__Rare_E.html
+tobit http://gking.harvard.edu/zelig/docs/_TT_tobit_TT__Regress.html
+tobit.bayes http://gking.harvard.edu/zelig/docs/_TT_tobit_bayes_TT.html
weibull http://gking.harvard.edu/zelig/docs/_TT_weibull_TT__Weibul.html
diff --git a/demo/00Index b/demo/00Index
index f172183..df0db6e 100644
--- a/demo/00Index
+++ b/demo/00Index
@@ -1,41 +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
-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
+beta Beta regression and simulation
+blogit Bivariate Logit regression and simulation
+bivariate.probit New model example: Bivariate probit MLE
+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
+normal.regression New model example for normal regression
+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
+tobit Regression for classical tobit 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/bivariate.probit.R b/demo/bivariate.probit.R
new file mode 100644
index 0000000..4a67865
--- /dev/null
+++ b/demo/bivariate.probit.R
@@ -0,0 +1,81 @@
+describe.bivariate.probit <- function() {
+ category <- "bivaraite.dichotomous"
+ package <- list(name = "mvtnorm",
+ version = "0.7")
+ mu <- list(equations = 2, # Systematic component
+ tagsAllowed = TRUE,
+ depVar = TRUE,
+ expVar = TRUE)
+ rho <- list(equations = 1, # Optional systematic component
+ tagsAllowed = FALSE, # Estimated as an ancillary
+ depVar = FALSE, # parameter by default
+ expVar = TRUE)
+ pars <- list(mu = mu, rho = rho)
+ list(category = category, parameters = pars)
+}
+
+zelig2bivariate.probit <- function(formula, model, data, M, ...) {
+ require(mvtnorm)
+ mf <- match.call(expand.dots = TRUE)
+ mf$formula <- parse.formula(formula, model)
+ mf$model <- mf$M <- NULL
+ mf[[1]] <- as.name("bivariate.probit")
+ as.call(mf)
+}
+
+bivariate.probit <- function(formula, data, start.val = NULL, ...) {
+
+ fml <- formula
+ D <- model.frame(fml, data = data)
+ X <- model.matrix(fml, data = D, eqn = c("mu1", "mu2")) # [1]
+ Xrho <- model.matrix(fml, data = D, eqn = "rho")
+
+ Y <- model.response(D)
+ terms <- attr(D,"terms")
+ start.val <- set.start(start.val, terms)
+ start.val <- put.start(start.val, 1, terms, eqn = "rho")
+
+ log.lik <- function(par, X, Y, terms) {
+ Beta <- parse.par(par, terms, eqn = c("mu1", "mu2")) # [2]
+ gamma <- parse.par(par, terms, eqn = "rho")
+ rho <- (exp(Xrho %*% gamma) - 1) / (1 + exp(Xrho %*% gamma))
+
+ mu <- X %*% Beta # [3]
+ llik <- 0
+ for (i in 1:nrow(mu)){
+ Sigma <- matrix(c(1, rho[i,], rho[i,], 1), 2, 2)
+ if (Y[i,1]==1)
+ if (Y[i,2]==1)
+ llik <- llik + log(pmvnorm(lower = c(0, 0), upper = c(Inf, Inf),
+ mean = mu[i,], corr = Sigma))
+ else
+ llik <- llik + log(pmvnorm(lower = c(0, -Inf), upper = c(Inf, 0),
+ mean = mu[i,], corr = Sigma))
+ else
+ if (Y[i,2]==1)
+ llik <- llik + log(pmvnorm(lower = c(-Inf, 0), upper = c(0, Inf),
+ mean = mu[i,], corr = Sigma))
+ else
+ llik <- llik + log(pmvnorm(lower = c(-Inf, -Inf), upper = c(0, 0),
+ mean = mu[i,], corr = Sigma))
+ }
+ return(llik)
+ }
+
+ res <- optim(start.val, log.lik, method = "BFGS",
+ hessian = TRUE, control = list(fnscale = -1),
+ X = X, Y = Y, terms = terms, ...)
+
+ fit <- model.end(res, D)
+ class(fit) <- "bivariate.probit"
+ fit
+}
+
+data(sanction)
+user.prompt()
+z.out1 <- zelig(cbind(import, export) ~ coop + cost + target, model = "bivariate.probit", data = sanction)
+user.prompt()
+z.out2 <- zelig(list(mu1 = import ~ coop, mu2 = export ~ cost + target),
+ model = "bivariate.probit", data = sanction)
+user.prompt()
+
diff --git a/demo/blogit.R b/demo/blogit.R
index dce4635..c541edd 100644
--- a/demo/blogit.R
+++ b/demo/blogit.R
@@ -38,8 +38,7 @@ user.prompt()
# Estimate the statistical model, with import a function of coop
# in the first equation and export a function of cost and target
# in the second equation, by using the zeros argument:
-z.out2 <- zelig(cbind(import, export) ~ coop + cost + target,
- omit = list("1" = "cost", "2" = c("coop", "target")),
+z.out2 <- zelig(list(mu1=import~coop,mu2=export~cost+target),
model = "blogit", data = sanction)
user.prompt()
print(summary(z.out2))
diff --git a/demo/bprobit.R b/demo/bprobit.R
index c5ee1c7..97d63b5 100644
--- a/demo/bprobit.R
+++ b/demo/bprobit.R
@@ -40,19 +40,18 @@ plot(s.out1)
# in the first equation and export a function of cost and target
# in the second equation, by using the zeros argument:
user.prompt()
-z.out2 <- zelig(cbind(import, export) ~ coop + cost + target,
- omit = list("1"="coop", "2"="cost", "3" = "target"),
+z.out2 <- zelig(list(mu1 = import ~ coop, mu2 = export ~ cost + target),
model = "bprobit", data = sanction)
user.prompt()
print(summary(z.out2))
# Set the explanatory variables to their default values:
user.prompt()
-Xval2 <- setx(z.out2)
+x.out2 <- setx(z.out2)
# Simulate draws from the posterior distribution:
user.prompt()
-s.out2 <- sim(z.out2, x = Xval2)
+s.out2 <- sim(z.out2, x = x.out2)
user.prompt()
print(summary(s.out2))
@@ -67,9 +66,8 @@ plot(s.out2)
# some or all explanatory variables, {\it and} the effects of the shared
# explanatory variables are jointly estimated. For example,
user.prompt()
-z.out3 <- zelig(cbind(import, export) ~ coop + cost + target,
- constrain = list("1"=c("coop", "cost", "target"),
- "2"=c("coop", "cost", "target")),
+z.out3 <- zelig(list(mu1 = import ~ tag(coop,"coop") + tag(cost, "cost") + tag (target, "target"),
+ mu2 = export ~ tag(coop,"coop") + tag(cost, "cost") + tag (target, "target")),
model = "bprobit", data = sanction)
user.prompt()
print(summary(z.out3))
diff --git a/demo/factor.mix.R b/demo/factor.mix.R
index 65964fd..45e4c70 100644
--- a/demo/factor.mix.R
+++ b/demo/factor.mix.R
@@ -4,7 +4,8 @@ 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)
+ burnin=5000,mcmc=100000, thin=50, verbose=TRUE,
+ L0=0.25,tune=1.2)
user.prompt()
## Checking for convergence before summarizing the estimates:
diff --git a/demo/factor.ord.R b/demo/factor.ord.R
index 0bc8321..21a9ea9 100644
--- a/demo/factor.ord.R
+++ b/demo/factor.ord.R
@@ -5,8 +5,8 @@ data(newpainters)
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)
+ burin=5000,mcmc=30000, thin=5, verbose=TRUE,
+ L0=0.5,tune=1.2)
user.prompt()
diff --git a/demo/irt1d.R b/demo/irt1d.R
index 09e898a..572c5e9 100644
--- a/demo/irt1d.R
+++ b/demo/irt1d.R
@@ -10,7 +10,7 @@ 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)
+ thin=20, verbose=FALSE)
user.prompt()
## Checking for convergence before summarizing the estimates:
diff --git a/demo/match.R b/demo/match.R
index 036e044..28676a2 100644
--- a/demo/match.R
+++ b/demo/match.R
@@ -105,6 +105,5 @@ summary(s.out3, subset = 2)
user.prompt()
summary(s.out3, subset = 3)
-user.prompt()
diff --git a/demo/mlogit.R b/demo/mlogit.R
index 91df41b..2328e88 100644
--- a/demo/mlogit.R
+++ b/demo/mlogit.R
@@ -1,27 +1,49 @@
data(mexico)
user.prompt()
-z.out <- zelig(vote88 ~ pristr + othcok + othsocok, model = "mlogit",
+z.out1 <- zelig(as.factor(vote88) ~ pristr + othcok + othsocok, model = "mlogit",
data = mexico)
user.prompt()
-print(summary(z.out))
+print(summary(z.out1))
user.prompt()
-x.weak <- setx(z.out, pristr = 1)
-x.strong <- setx(z.out, pristr = 3)
+x.weak <- setx(z.out1, pristr = 1)
+x.strong <- setx(z.out1, pristr = 3)
user.prompt()
-s.out <- sim(z.out, x = x.strong, x1 = x.weak)
+s.out1 <- sim(z.out1, x = x.strong, x1 = x.weak)
user.prompt()
-print(summary(s.out))
+print(summary(s.out1))
user.prompt()
-ev.weak <- s.out$qi$ev + s.out$qi$fd
+ev.weak <- s.out1$qi$ev + s.out1$qi$fd
user.prompt()
-library(vcd)
-ternaryplot(s.out$qi$ev, pch = ".", col = "blue",
+ternaryplot(s.out1$qi$ev, pch = ".", col = "blue",
main = "1988 Mexican Presidential Election")
user.prompt()
ternarypoints(ev.weak, pch = ".", col = "red")
+# Specifying different sets of explanatory variables for each factor level
+user.prompt()
+z.out2 <- zelig(list(id(vote88,"1")~pristr + othcok, id(vote88,"2")~othsocok), model = "mlogit",
+ data = mexico)
+user.prompt()
+print(summary(z.out2))
+
+user.prompt()
+x.weak <- setx(z.out2, pristr = 1)
+x.strong <- setx(z.out2, pristr = 3)
+
+user.prompt()
+s.out2 <- sim(z.out2, x = x.strong, x1 = x.weak)
+user.prompt()
+print(summary(s.out2))
+
+user.prompt()
+ev.weak <- s.out2$qi$ev + s.out2$qi$fd
+user.prompt()
+ternaryplot(s.out2$qi$ev, pch = ".", col = "blue",
+ main = "1988 Mexican Presidential Election")
+user.prompt()
+ternarypoints(ev.weak, pch = ".", col = "red")
diff --git a/demo/normal.regression.R b/demo/normal.regression.R
new file mode 100644
index 0000000..b35ddd6
--- /dev/null
+++ b/demo/normal.regression.R
@@ -0,0 +1,102 @@
+describe.normal.regression <- function() {
+ category <- "continuous"
+ mu <- list(equations = 1, # Systematic component
+ tagsAllowed = FALSE,
+ depVar = TRUE,
+ expVar = TRUE)
+ sigma2 <- list(equations = 1, # Scalar ancillary parameter
+ tagsAllowed = FALSE,
+ depVar = FALSE,
+ expVar = FALSE)
+ pars <- list(mu = mu, sigma2 = sigma2)
+ model <- list(category = category, parameters = pars)
+}
+
+zelig2normal.regression <- function(formula, model, data, M, ...) {
+ mf <- match.call(expand.dots = TRUE) # [1]
+ mf$model <- mf$M <- NULL # [2]
+ mf[[1]] <- as.name("normal.regression") # [3]
+ as.call(mf) # [4]
+}
+
+normal.regression <- function(formula, data, start.val = NULL, ...) {
+
+ fml <- parse.formula(formula, model = "normal.regression") # [1]
+ D <- model.frame(fml, data = data)
+ X <- model.matrix(fml, data = D)
+ Y <- model.response(D)
+ terms <- attr(D, "terms")
+ n <- nrow(X)
+
+ start.val <- set.start(start.val, terms)
+ start.val <- put.start(start.val, 1, terms, eqn = "sigma2") # [2]
+
+ ll.normal <- function(par, X, Y, n, terms) { # [3]
+ beta <- parse.par(par, terms, eqn = "mu")
+ sigma2 <- exp(parse.par(par, terms, eqn = "sigma2"))
+ -n / 2 * log(2 * pi * sigma2) -
+ 1/(2 * sigma2) * sum((Y - X %*% beta)^2)
+ }
+
+ res <- optim(start.val, ll.normal, method = "BFGS",
+ hessian = TRUE, control = list(fnscale = -1),
+ X = X, Y = Y, n = n, terms = terms, ...) # [4]
+
+ fit <- model.end(res, D) # [5]
+ fit$n <- n
+ class(fit) <- "normal" # [6]
+ fit # [7]
+}
+
+param.normal <- function(object, num = NULL, bootstrap = FALSE,
+ terms = NULL) {
+ if (!bootstrap) {
+ par <- mvrnorm(num, mu = coef(object), Sigma = vcov(object))
+ Beta <- parse.par(par, terms = terms, eqn = "mu")
+ sigma2 <- exp(parse.par(par, terms = terms, eqn = "sigma2"))
+ res <- cbind(Beta, sigma2)
+ }
+ else {
+ par <- coef(object)
+ Beta <- parse.par(par, terms = terms, eqn = "mu")
+ sigma2 <- exp(parse.par(par, terms = terms, eqn = "sigma2"))
+ res <- c(coef, sigma2)
+ }
+ res
+}
+
+qi.normal <- function(object, par, x, x1 = NULL, y = NULL) {
+ Beta <- parse.par(par, eqn = "mu")
+ sigma2 <- parse.par(par, eqn = "sigma2")
+ ev <- Beta %*% t(x)
+ pr <- matrix(NA, ncol = ncol(ev), nrow = nrow(ev))
+ for (i in 1:ncol(ev)) # Using R's built-in poisson generator.
+ pr[,i] <- rnorm(length(ev[,i]), mean = ev[,i], sigma = sd(sigma2[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)){
+ ev1 <- par %*% t(x1)
+ qi$fd <- ev1 - ev
+ qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
+ }
+ if (!is.null(y)) {
+ yvar <- matrix(rep(y, nrow(par)), nrow = nrow(par), byrow = TRUE)
+ tmp.ev <- yvar - qi$ev
+ tmp.pr <- yvar - qi$pr
+ qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(par))
+ qi$ate.pr <- matrix(apply(tmp.pr, 1, mean), nrow = nrow(par))
+ 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)
+}
+
+data(macro)
+
+user.prompt()
+
+z.out <- zelig(unem ~ gdp + capmob + trade, model = "normal.regression",
+data = macro)
+x.out <- setx(z.out)
+s.out <- setx(z.out, x = x.out)
diff --git a/demo/oprobit.bayes.R b/demo/oprobit.bayes.R
index 7682503..b1a5a03 100644
--- a/demo/oprobit.bayes.R
+++ b/demo/oprobit.bayes.R
@@ -9,7 +9,7 @@ sanction$ncost <- factor(sanction$ncost, ordered = TRUE,
## Estimating the model using oprobit.bayes:
z.out <- zelig(ncost ~ mil + coop, model = "oprobit.bayes",
- data = sanction, verbose=TRUE)
+ data = sanction, verbose=FALSE, tune=0.3)
user.prompt()
diff --git a/demo/tobit.R b/demo/tobit.R
new file mode 100644
index 0000000..5a8ec98
--- /dev/null
+++ b/demo/tobit.R
@@ -0,0 +1,47 @@
+## Attaching the example dataset:
+data(tobin)
+
+## Estimating the model using tobit.bayes:
+z.out <- zelig(durable ~ age + quant, data = tobin, model = "tobit")
+user.prompt()
+
+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/current.packages.Rd b/man/current.packages.Rd
new file mode 100644
index 0000000..d8a65a3
--- /dev/null
+++ b/man/current.packages.Rd
@@ -0,0 +1,36 @@
+\name{current.packages}
+
+\alias{current.packages}
+
+\title{Find all packages in a dependency chain}
+
+\description{
+ Use \code{current.packages} to find all the packages suggested or
+required by a given package, and the currently installed version number
+for each.
+}
+
+\usage{
+current.packages(package)
+}
+
+\arguments{
+\item{package}{a character string corresponding to the name of an
+installed package}
+}
+
+\value{
+ A matrix containing the current version number of the packages
+suggested or required by \code{package}.
+}
+
+\seealso{\code{packageDescription}}
+
+\author{ Olivia Lau <\email{olau at fas.harvard.edu}>
+}
+
+\examples{
+current.packages("Zelig")
+}
+
+\keyword{file}
diff --git a/man/model.end.Rd b/man/model.end.Rd
new file mode 100644
index 0000000..c004d88
--- /dev/null
+++ b/man/model.end.Rd
@@ -0,0 +1,41 @@
+\name{model.end}
+\alias{model.end}
+\title{Cleaning up after optimization}
+
+\description{ The \code{model.end} function creates a list of regression output from \code{\link{optim}}
+output. The list includes coefficients (from the \code{\link{optim}} \code{par} output), a
+variance-covariance matrix (from the \code{\link{optim}} Hessian output), and any terms, contrasts, or
+xlevels (from the model frame). Use \code{model.end} after calling \code{\link{optim}}, but before
+assigning a
+class to the regression output.}
+
+\usage{
+model.end(res, mf)
+}
+
+\arguments{
+\item{res}{the output from \code{\link{optim}} or another fitting-algorithm}
+\item{mf}{the model frame output by \code{\link{model.frame}}}
+}
+
+\value{
+A list of regression output, including:
+\itemize{
+\item{coefficients}{the optimized parameters}
+\item{variance}{the variance-covariance matrix (the negative
+ inverse of the Hessian matrix returned from the optimization
+ procedure)}
+\item{terms}{the terms object. See \code{\link{terms.object}}
+ for more information}
+\item{\dots}{additional elements passed from \code{res}}
+}}
+
+\seealso{The full Zelig manual at \url{http://gking.harvard.edu/zelig} for examples.}
+
+\author{
+ Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Olivia Lau <\email{olau at fas.harvard.edu}>; Ferdinand Alimadhi
+<\email{falimadhi at iq.harvard.edu}>
+}
+
+\keyword{utilities}
diff --git a/man/model.frame.multiple.Rd b/man/model.frame.multiple.Rd
new file mode 100644
index 0000000..9075a43
--- /dev/null
+++ b/man/model.frame.multiple.Rd
@@ -0,0 +1,60 @@
+\name{model.frame.multiple}
+\alias{model.frame.multiple}
+\title{Extracting the ``environment'' of a model formula}
+
+\description{ Use \code{model.frame.multiple} after \code{\link{parse.par}} to create a
+data frame of the unique variables identified in the formula (or list
+of formulas).}
+
+\usage{
+model.frame.multiple(formula, data, ...)
+}
+
+\arguments{
+\item{formula}{a list of formulas of class \code{"multiple"}, returned from \code{\link{parse.par}}}
+\item{data}{a data frame containing all the variables used in \code{formula}}
+\item{\dots}{additional arguments passed to \code{\link{model.frame.default}}}
+}
+
+\value{
+The output is a data frame (with a terms attribute) containing all the
+unique explanatory and response variables identified in the list of
+formulas. By default, missing (\code{NA}) values are listwise deleted.
+
+If \code{as.factor} appears on the left-hand side, the response
+variables will be returned as an indicator (0/1) matrix with columns
+corresponding to the unique levels in the factor variable.
+
+If any formula contains more than one \code{tag} statement, \code{model.frame.multiple}
+will return the original variable in the data frame and use the \code{tag} information in the terms
+attribute only.
+}
+
+\examples{
+\dontrun{
+data(sanction)
+formulae <- list(import ~ coop + cost + target,
+ export ~ coop + cost + target)
+fml <- parse.formula(formulae, model = "bivariate.logit")
+D <- model.frame(fml, data = sanction)
+}}
+
+\seealso{\code{\link{model.matrix.default}}, \code{\link{parse.formula}} and the full Zelig manual at
+ \url{http://gking.harvard.edu/zelig}}
+
+\author{
+ Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Olivia Lau <\email{olau at fas.harvard.edu}>; Ferdinand Alimadhi
+<\email{falimadhi at iq.harvard.edu}>
+}
+
+\keyword{utilities}
+
+
+
+
+
+
+
+
+
diff --git a/man/model.matrix.multiple.Rd b/man/model.matrix.multiple.Rd
new file mode 100644
index 0000000..9d595e0
--- /dev/null
+++ b/man/model.matrix.multiple.Rd
@@ -0,0 +1,91 @@
+\name{model.matrix.multiple}
+\alias{model.matrix.multiple}
+\title{Design matrix for multivariate models}
+
+\description{Use \code{model.matrix.multiple} after \code{\link{parse.formula}} to
+create a design matrix for multiple-equation models. }
+
+\usage{
+model.matrix.multiple(object, data, shape = "compact", eqn = NULL, ...)
+}
+
+\arguments{
+\item{object}{the list of formulas output from \code{\link{parse.formula}}}
+\item{data}{a data frame created with \code{\link{model.frame.multiple}}}
+\item{shape}{a character string specifying the shape of the outputed matrix. Available options are
+\itemize{
+\item{"compact"}{(default) the output matrix will be an \eqn{n \times v}{n x v},
+where \eqn{v} is the number of unique variables in all of the equations
+(including the intercept term)}
+\item{"array"}{the output is an \eqn{n \times K \times J}{n x K x J} array where \eqn{J} is the
+total number of equations and \eqn{K} is the total number of parameters
+across all the equations. If a variable is not in a certain equation,
+it is observed as a vector of 0s.}
+\item{"stacked"}{the output will be a \eqn{2n \times K}{2n x K} matrix where \eqn{K} is the total number of
+parameters across all the equations.}
+}}
+\item{eqn}{a character string or a vector of character strings identifying the equations from which to
+construct the design matrix. The defaults to \code{NULL}, which only uses the systematic
+ parameters (for which \code{DepVar = TRUE} in the appropriate \code{describe.model} function)}
+\item{\dots}{additional arguments passed to \code{model.matrix.default}}
+}
+
+\value{A design matrix or array, depending on the options chosen in \code{shape}, with appropriate terms
+attributes.}
+
+\examples{
+
+# Let's say that the name of the model is "bivariate.probit", and
+# the corresponding describe function is describe.bivariate.probit(),
+# which identifies mu1 and mu2 as systematic components, and an
+# ancillary parameter rho, which may be parameterized, but is estimated
+# as a scalar by default. Let par be the parameter vector (including
+# parameters for rho), formulae a user-specified formula, and mydata
+# the user specified data frame.
+
+# Acceptable combinations of parse.par() and model.matrix() are as follows:
+## Setting up
+\dontrun{
+data(sanction)
+formulae <- cbind(import, export) ~ coop + cost + target
+fml <- parse.formula(formulae, model = "bivariate.probit")
+D <- model.frame(fml, data = sanction)
+terms <- attr(D, "terms")
+
+## Intuitive option
+Beta <- parse.par(par, terms, shape = "vector", eqn = c("mu1", "mu2"))
+X <- model.matrix(fml, data = D, shape = "stacked", eqn = c("mu1", "mu2")
+eta <- X %*% Beta
+
+## Memory-efficient (compact) option (default)
+Beta <- parse.par(par, terms, eqn = c("mu1", "mu2"))
+X <- model.matrix(fml, data = D, eqn = c("mu1", "mu2"))
+eta <- X %*% Beta
+
+## Computationally-efficient (array) option
+Beta <- parse.par(par, terms, shape = "vector", eqn = c("mu1", "mu2"))
+X <- model.matrix(fml, data = D, shape = "array", eqn = c("mu1", "mu2"))
+eta <- apply(X, 3, '%*%', Beta)
+}}
+
+\seealso{\code{\link{parse.par}}, \code{\link{parse.formula}} and the full Zelig manual at
+ \url{http://gking.harvard.edu/zelig}}
+
+\author{
+ Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Olivia Lau <\email{olau at fas.harvard.edu}>; Ferdinand Alimadhi
+<\email{falimadhi at iq.harvard.edu}>
+}
+
+\keyword{utilities}
+
+
+
+
+
+
+
+
+
+
+
diff --git a/man/parse.formula.Rd b/man/parse.formula.Rd
new file mode 100644
index 0000000..e6b7354
--- /dev/null
+++ b/man/parse.formula.Rd
@@ -0,0 +1,88 @@
+\name{parse.formula}
+\alias{parse.formula}
+
+\title{Parsing user-input formulas into multiple syntax}
+
+\description{Parse the input formula (or list of formulas) into the
+standard format described below. Since labels for this format will vary
+by model, \code{parse.formula} will evaluate a function \code{describe.model},
+where \code{model} is given as an input to \code{parse.formula}.
+
+If the \code{describe.model} function has more than one parameter for
+which \code{ExpVar = TRUE} and \code{DepVar = TRUE}, then the
+user-specified equations must have labels to match those parameters,
+else \code{parse.formula} should return an error. In addition, if the
+formula entries are not unambiguous, then \code{parse.formula} returns an error.
+}
+
+\usage{
+parse.formula(formula, model, data = NULL)
+}
+
+\arguments{
+\item{formula}{either a single formula or a list of \code{formula} objects}
+\item{model}{a character string specifying the name of the model}
+\item{data}{an optional data frame for models that require a factor response variable}
+}
+
+\value{The output is a list of formula objects with class
+\code{c("multiple", "list")}. Let's say that the name of the model is
+\code{"bivariate.probit"}, and the corresponding describe function is
+\code{describe.bivariate.probit}, which identifies \code{mu1} and
+\code{mu2} as systematic components, and an ancillary parameter \code{rho}, which
+may be parameterized, but is estimated as a scalar by default.
+}
+
+\details{Acceptable user inputs are as follows:
+
+\tabular{lll}{
+ \tab User Input \tab Output from \code{parse.formula}\cr
+\tab \tab \cr
+Same covariates, \tab cbind(y1, y2) ~ x1 + x2 * x3 \tab list(mu1 = y1 ~ x1 + x2 * x3,\cr
+separate effects \tab \tab mu2 = y2 ~ x1 + x2 * x3,\cr
+ \tab \tab rho = ~ 1)\cr
+\tab \tab \cr
+With \code{rho} as a \tab list(cbind(y1, y2) ~ x1 + x2, \tab list(mu1 = y1 ~ x1 + x2,\cr
+systematic equation \tab rho = ~ x4 + x5) \tab mu2 = y2 ~ x1 + x2,\cr
+ \tab \tab rho = ~ x4 + x5)\cr
+\tab \tab \cr
+With constraints \tab list(mu1 = y1 ~ x1 + tag(x2, "x2"), \tab list(mu1 = y1 ~ x1 + tag(x2, "x2"),\cr
+(same variable) \tab mu2 = y2 ~ x3 + tag(x2, "x2")) \tab mu2 = y2 ~ x3 + tag(x2, "x2"),\cr
+ \tab \tab rho = ~ 1)\cr
+\tab \tab \cr
+With constraints \tab list(mu1 = y1 ~ x1 + tag(x2, "z1"), \tab list(mu1 = y1 ~ x1 + tag(x2, "z1"),\cr
+(different variables) \tab mu2 = y2 ~ x3 + tag(x4, "z1")) \tab mu2 = y2 ~ x3 + tag(x4, "z1"),\cr
+ \tab \tab rho = ~ 1)\cr
+}}
+
+\examples{
+\dontrun{
+data(sanction)
+formulae <- list(cbind(import, export) ~ coop + cost + target)
+fml <- parse.formula(formulae, model = "bivariate.probit")
+D <- model.frame(fml, data = sanction)
+}}
+
+\seealso{
+\code{\link{parse.par}}, \code{\link{model.frame.multiple}},
+\code{\link{model.matrix.multiple}}, and the full Zelig manual at
+ \url{http://gking.harvard.edu/zelig}.
+}
+
+\author{
+ Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Olivia Lau <\email{olau at fas.harvard.edu}>; Ferdinand Alimadhi
+<\email{falimadhi at iq.harvard.edu}>
+}
+
+\keyword{utilities}
+
+
+
+
+
+
+
+
+
+
diff --git a/man/parse.par.Rd b/man/parse.par.Rd
new file mode 100644
index 0000000..5d97a52
--- /dev/null
+++ b/man/parse.par.Rd
@@ -0,0 +1,82 @@
+\name{parse.par}
+\alias{parse.par}
+
+\title{Select and reshape parameter vectors}
+
+\description{ The \code{parse.par} function reshapes parameter vectors for
+comfortability with the output matrix from \code{\link{model.matrix.multiple}}.
+Use \code{parse.par} to identify sets of parameters; for example, within
+optimization functions that require vector input, or within \code{qi}
+functions that take matrix input of all parameters as a lump.
+}
+
+\usage{
+parse.par(par, terms, shape = "matrix", eqn = NULL)
+}
+
+\arguments{
+\item{par}{the vector (or matrix) of parameters}
+\item{terms}{the terms from either \code{\link{model.frame.multiple}} or
+\code{\link{model.matrix.multiple}}}
+\item{shape}{a character string (either \code{"matrix"} or \code{"vector"})
+that identifies the type of output structure}
+\item{eqn}{a character string (or strings) that identify the
+parameters that you would like to subset from the larger \code{par}
+structure}
+}
+
+\value{
+A matrix or vector of the sub-setted (and reshaped) parameters for the specified
+parameters given in \code{"eqn"}. By default, \code{eqn = NULL}, such that all systematic
+components are selected. (Systematic components have \code{ExpVar = TRUE} in the appropriate
+\code{describe.model} function.)
+
+If an ancillary parameter (for which \code{ExpVar = FALSE} in
+\code{describe.model}) is specified in \code{eqn}, it is
+always returned as a vector (ignoring \code{shape}). (Ancillary
+parameters are all parameters that have intercept only formulas.)
+}
+\examples{
+# Let's say that the name of the model is "bivariate.probit", and
+# the corresponding describe function is describe.bivariate.probit(),
+# which identifies mu1 and mu2 as systematic components, and an
+# ancillary parameter rho, which may be parameterized, but is estimated
+# as a scalar by default. Let par be the parameter vector (including
+# parameters for rho), formulae a user-specified formula, and mydata
+# the user specified data frame.
+
+# Acceptable combinations of parse.par() and model.matrix() are as follows:
+## Setting up
+\dontrun{
+data(sanction)
+formulae <- cbind(import, export) ~ coop + cost + target
+fml <- parse.formula(formulae, model = "bivariate.probit")
+D <- model.frame(fml, data = sanction)
+terms <- attr(D, "terms")
+
+## Intuitive option
+Beta <- parse.par(par, terms, shape = "vector", eqn = c("mu1", "mu2"))
+X <- model.matrix(fml, data = D, shape = "stacked", eqn = c("mu1", "mu2")
+eta <- X %*% Beta
+
+## Memory-efficient (compact) option (default)
+Beta <- parse.par(par, terms, eqn = c("mu1", "mu2"))
+X <- model.matrix(fml, data = D, eqn = c("mu1", "mu2"))
+eta <- X %*% Beta
+
+## Computationally-efficient (array) option
+Beta <- parse.par(par, terms, shape = "vector", eqn = c("mu1", "mu2"))
+X <- model.matrix(fml, data = D, shape = "array", eqn = c("mu1", "mu2"))
+eta <- apply(X, 3, '%*%', Beta)
+}}
+
+\seealso{\code{\link{model.matrix.multiple}}, \code{\link{parse.formula}} and the full Zelig manual at
+ \url{http://gking.harvard.edu/zelig}}
+
+\author{
+ Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Olivia Lau <\email{olau at fas.harvard.edu}>; Ferdinand Alimadhi
+<\email{falimadhi at iq.harvard.edu}>
+}
+
+\keyword{utilities}
diff --git a/man/plot.ci.Rd b/man/plot.ci.Rd
index e855b0c..204af5e 100644
--- a/man/plot.ci.Rd
+++ b/man/plot.ci.Rd
@@ -41,7 +41,7 @@ variable. You may save this plot using the commands described in the
Zelig manual (\url{http://gking.harvard.edu/zelig}).
}
-\example{
+\examples{
data(turnout)
z.out <- zelig(vote ~ race + educate + age + I(age^2) + income,
model = "logit", data = turnout)
diff --git a/man/put.start.Rd b/man/put.start.Rd
new file mode 100644
index 0000000..2a41bd2
--- /dev/null
+++ b/man/put.start.Rd
@@ -0,0 +1,34 @@
+\name{put.start}
+\alias{put.start}
+\title{Set specific starting values for certain parameters}
+
+\description{ After calling \code{\link{set.start}} to create default starting values, use \code{put.start}
+to change starting values for specific parameters or parameter sets. }
+
+\usage{
+put.start(start.val, value, terms, eqn)
+}
+
+\arguments{
+\item{start.val}{the vector of starting values created by \code{\link{set.start}}}
+\item{value}{the scalar or vector of replacement starting values}
+\item{terms}{the terms output from \code{\link{model.frame.multiple}}}
+\item{eqn}{character vector of the parameters for which you would like to replace
+the default values with \code{value}}
+}
+
+\value{A vector of starting values (of the same length as \code{start.val})}
+
+\seealso{\code{\link{set.start}}, and the full Zelig manual at
+ \url{http://gking.harvard.edu/zelig}.
+}
+
+\author{
+ Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Olivia Lau <\email{olau at fas.harvard.edu}>; Ferdinand Alimadhi
+<\email{falimadhi at iq.harvard.edu}>
+}
+
+
+\keyword{utilities}
+
diff --git a/man/repl.Rd b/man/repl.Rd
index f97fd0c..72e16e2 100644
--- a/man/repl.Rd
+++ b/man/repl.Rd
@@ -68,9 +68,9 @@ this error decreases.
\code{\link{sim}}. In addition, the full Zelig manual may be
accessed online at \url{http://gking.harvard.edu/zelig}. }
-\example{
+\examples{
data(turnout)
-z.out <- zelig(vote ~ race + educate, model = "logit", data = turnout[1:1000,$
+z.out <- zelig(vote ~ race + educate, model = "logit", data = turnout[1:1000,])
x.out <- setx(z.out)
s.out <- sim(z.out, x = x.out)
z.rep <- repl(z.out)
diff --git a/man/set.start.Rd b/man/set.start.Rd
new file mode 100644
index 0000000..64e8ae6
--- /dev/null
+++ b/man/set.start.Rd
@@ -0,0 +1,40 @@
+\name{set.start}
+\alias{set.start}
+\title{Set starting values for all parameters}
+
+\description{After using \code{\link{parse.par}} and \code{\link{model.matrix.multiple}}, use
+\code{set.start} to set starting values for all parameters. By default, starting values are set to 0. If
+you wish to select alternative starting values for certain parameters, use \code{\link{put.start}} after
+\code{set.start}.}
+
+\usage{
+set.start(start.val = NULL, terms)
+}
+
+\arguments{
+\item{start.val}{user-specified starting values. If \code{NULL} (default), the default
+starting values for all parameters are set to 0.}
+\item{terms}{the terms output from \code{\link{model.frame.multiple}}}
+}
+
+\value{
+A named vector of starting values for all parameters specified in \code{terms}, defaulting to 0.
+}
+
+\examples{
+\dontrun{
+fml <- parse.formula(formula, model = "bivariate.probit")
+D <- model.frame(fml, data = data)
+terms <- attr(D, "terms")
+start.val <- set.start(start.val = NULL, terms)
+}}
+
+\seealso{\code{\link{put.start}}, \code{\link{parse.par}}, \code{\link{model.frame.multiple}}, and the
+full Zelig manual at \url{http://gking.harvard.edu/zelig}.}
+
+\author{
+ Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
+ <\email{king at harvard.edu}>; Olivia Lau <\email{olau at fas.harvard.edu}>; Ferdinand Alimadhi
+<\email{falimadhi at iq.harvard.edu}>
+}
+\keyword{utilities}
diff --git a/man/setx.Rd b/man/setx.Rd
index 9860462..b4c0e22 100644
--- a/man/setx.Rd
+++ b/man/setx.Rd
@@ -57,9 +57,7 @@ x.out <- setx(object, fn = list(numeric = mean, ordered = median,
\code{setx} returns the selected values calculated over the entire
data frame. If you wish to calculate values over just one subset of
the data frame, the 5th subset for example, you may use:
-\begin{verbatim}
-x.out <- setx(z.out[[5]])
-\end{verbatim}
+\code{x.out <- setx(z.out[[5]])}
For conditional prediction, \code{x.out} includes the model matrix
and the dependent variables. For multiple analyses (when choosing
@@ -67,7 +65,7 @@ For conditional prediction, \code{x.out} includes the model matrix
observed explanatory variables in each subset.
}
-\example{
+\examples{
# Unconditional prediction:
data(turnout)
z.out <- zelig(vote ~ race + educate, model = "logit", data = turnout)
@@ -85,20 +83,24 @@ x.out <- setx(z.out, data = turnout[1001:2000,])
s.out <- sim(z.out, x = x.out)
# Using a user-defined function in fn:
+\dontrun{
quants <- function(x)
quantile(x, 0.25)
x.out <- setx(z.out, fn = list(numeric = quants))
+}
# Conditional prediction:
+\dontrun{library(MatchIt)
data(lalonde)
match.out <- matchit(treat ~ age + educ + black + hispan + married +
nodegree + re74 + re75, data = lalonde)
-z.out <- zelig(re78 ~ pscore, data = match.data(match, "control"),
+z.out <- zelig(re78 ~ distance, data = match.data(match.out, "control"),
model = "ls")
x.out <- setx(z.out, fn = NULL, data = match.data(match.out, "treat"),
cond = TRUE)
s.out <- sim(z.out, x = x.out)
}
+}
\seealso{ The full Zelig manual may be accessed online at
\url{http://gking.harvard.edu/zelig}. }
diff --git a/man/ternaryplot.Rd b/man/ternaryplot.Rd
new file mode 100644
index 0000000..74275a3
--- /dev/null
+++ b/man/ternaryplot.Rd
@@ -0,0 +1,72 @@
+\name{ternaryplot}
+\alias{ternaryplot}
+%- Also NEED an `\alias' for EACH other topic documented here.
+\title{Ternary diagram}
+\description{
+Visualizes compositional, 3-dimensional data in an equilateral triangle
+(from the vcd library, Version 0.1-3.3, Date 2004-04-21), using plot graphics.
+Differs from implementation in vcd (0.9-7), which uses grid graphics.}
+\usage{
+ternaryplot(x, scale = 1, dimnames = NULL, dimnames.position = c("corner","edge","none"),
+ dimnames.color = "black", id = NULL, id.color = "black", coordinates = FALSE,
+ grid = TRUE, grid.color = "gray", labels = c("inside", "outside", "none"),
+ labels.color = "darkgray", border = "black", bg = "white", pch = 19, cex = 1,
+ prop.size = FALSE, col = "red", main = "ternary plot", ...)
+}
+\arguments{
+ \item{x}{a matrix with three columns.}
+ \item{scale}{row sums scale to be used.}
+ \item{dimnames}{dimension labels (defaults to the column names of
+ \code{x}).}
+ \item{dimnames.position, dimnames.color}{position and color of dimension labels.}
+ \item{id}{optional labels to be plotted below the plot
+ symbols. \code{coordinates} and \code{id} are mutual exclusive.}
+ \item{id.color}{color of these labels.}
+ \item{coordinates}{if \code{TRUE}, the coordinates of the points are
+ plotted below them. \code{coordinates} and \code{id} are mutual exclusive.}
+ \item{grid}{if \code{TRUE}, a grid is plotted. May optionally
+ be a string indicating the line type (default: \code{"dotted"}).}
+ \item{grid.color}{grid color.}
+ \item{labels, labels.color}{position and color of the grid labels.}
+ \item{border}{color of the triangle border.}
+ \item{bg}{triangle background.}
+ \item{pch}{plotting character. Defaults to filled dots.}
+ \item{cex}{a numerical value giving the amount by which plotting text
+ and symbols should be scaled relative to the default. Ignored for
+ the symbol size if \code{prop.size} is not \code{FALSE}.}
+ \item{prop.size}{if \code{TRUE}, the symbol size is plotted
+ proportional to the row sum of the three variables, i.e. represents
+ the weight of the observation.}
+ \item{col}{plotting color.}
+ \item{main}{main title.}
+ \item{\dots}{additional graphics parameters (see \code{par})}
+}
+\details{
+A points' coordinates are found by computing the gravity center
+of mass points using the data entries as weights. Thus, the coordinates
+of a point P(a,b,c), \eqn{a + b + c = 1}, are: P(b + c/2, c * sqrt(3)/2).
+}
+
+\examples{
+data(mexico)
+z.out <- zelig(as.factor(vote88) ~ pristr + othcok + othsocok,
+ model = "mlogit", data = mexico)
+x.out <- setx(z.out)
+s.out <- sim(z.out, x = x.out)
+
+ternaryplot(s.out$qi$ev, pch = ".", col = "blue",
+ main = "1988 Mexican Presidential Election")
+}
+
+\seealso{\code{\link{ternarypoints}}}
+
+\references{
+M. Friendly (2000),
+\emph{Visualizing Categorical Data}. SAS Institute, Cary, NC.
+}
+\author{
+ David Meyer\cr
+ \email{david.meyer at ci.tuwien.ac.at}
+}
+
+\keyword{hplot}
diff --git a/man/ternarypoints.Rd b/man/ternarypoints.Rd
index 0488891..2e5cd60 100644
--- a/man/ternarypoints.Rd
+++ b/man/ternarypoints.Rd
@@ -31,13 +31,9 @@ ternary diagram. Use \code{ternaryplot} in the \code{vcd} library to
generate the main ternary diagram.
}
-\example{
-
-}
-
\seealso{ The full Zelig manual at
\url{http://gking.harvard.edu/zelig}, \code{points}, and
- \code{ternaryplot} in the \code{vcd} library. }
+ \code{\link{ternaryplot}}. }
\author{
Kosuke Imai <\email{kimai at princeton.edu}>; Gary King
diff --git a/man/zeligVDC.Rd b/man/zeligVDC.Rd
new file mode 100644
index 0000000..2642187
--- /dev/null
+++ b/man/zeligVDC.Rd
@@ -0,0 +1,82 @@
+\name{zeligDescribeModelXML}
+\alias{zeligDescribeModelXML}
+\alias{zeligInstalledModels}
+\alias{zeligListModels}
+\alias{zeligModelDependency}
+\alias{zeligGetSpecial}
+\title{ Zelig interface functions}
+\description{
+ Zelig interface functions. Used by VDC DSB to communicate with Zelig.
+}
+\usage{
+ zeligDescribeModelXML(modelName,force=F,schemaVersion="1.1")
+ zeligInstalledModels(inZeligOnly=T,schemaVersion="1.1")
+ zeligListModels(inZeligOnly=T)
+ zeligModelDependency(modelName,repos)
+ zeligGetSpecial(modelName)
+}
+
+\arguments{
+ \item{modelName}{Name of model as returned by zeligInstalledModels or zeligListModels.}
+ \item{inZeligOnly}{Flag, include only models in official Zelig distribution}
+ \item{repos}{URL of default repository to use}
+ \item{schemaVersion}{version of Zelig schema}
+ \item{force}{generate a description even if no custom description supplied}
+}
+
+\value{
+Use zeligInstalledModels and zeligListModels to determine what models are available in zelig
+for a particular schema level. Use zmodel2string(zeligDescribeModel()) to generate an XML
+instance describing a model. Use zeligModelDependencies to generate a list of package
+dependencies for models. Use zeligGetSpecial to get the name special function, if any,
+to apply to the outcome variables. All functions return NULL if results are
+not available for that model.
+}
+
+\examples{
+ # show all available models
+ zeligListModels(inZeligOnly=FALSE)
+ # show installed models
+ zeligInstalledModels()
+ # show dependency for normal.bayes
+ zeligModelDependency("normal.bayes","http://cran.r-project.org/")
+ # description of logit
+ cat(zeligDescribeModelXML("ologit"))
+ # special function for factor analysis
+ zeligGetSpecial("factor.mix")
+
+\dontshow{
+
+# test model lists
+zd= zeligInstalledModels(schemaVersion="1.1")
+if (length(zd)<8 || sum(zd=="ls")!=1 || length(zeligListModels())<25 ) {
+ stop("Failed zeligListModels/zeligInstalledModels self test")
+}
+
+if (zeligModelDependency("poisson.bayes","")[1]!="MCMCpack") {
+ stop("Failed zeligModelDependency self test")
+}
+
+if (zeligGetSpecial("factor.mix")!="cbind") {
+ stop("Failed zeligGetSpecial self test")
+}
+
+if (grep("explanatory", zeligDescribeModelXML("ologit"))!=1) {
+ stop("Failed zmodel2string/zeligDescribeModel self test")
+}
+
+}
+
+}
+
+\author{
+Micah Altman
+\email{thedata-users\@lists.sourceforge.net}
+\url{http://thedata.org}
+}
+
+
+\seealso{ \link[pkg:Zelig]{zelig}}
+
+\keyword{IO}
+\keyword{print}
--
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