[r-cran-zelig] 14/102: Import Upstream version 2.4-4
Andreas Tille
tille at debian.org
Sun Jan 8 16:58:09 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-zelig.
commit 840ea22c6de4d2103ebfb9abee2d3ff5761a9efd
Author: Andreas Tille <tille at debian.org>
Date: Sun Jan 8 09:39:00 2017 +0100
Import Upstream version 2.4-4
---
DESCRIPTION | 6 +-
NAMESPACE | 2 +-
R/MULTIPLE_NEW/make.parameters.R | 39 ++++
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/put.start.R | 18 ++
R/MULTIPLE_NEW/set.start.R | 15 ++
R/MULTIPLE_NEW/tag.R | 5 +
R/MULTIPLE_NEW/terms.list.R | 113 ++++++++++++
R/MULTIPLE_NEW/toBuildFormula.R | 17 ++
R/generate.start.R | 9 -
R/parse.par.R | 7 -
R/print.summary.MCMCZelig.R | 2 +-
R/print.summary.zelig.R | 2 +-
R/qi.MCMCZelig.R | 368 +++++++++++++------------------------
R/qi.survreg.R | 53 +++---
R/setx.default.R | 2 +-
README | 8 +
data/MatchIt.url.tab | 11 +-
data/Zelig.url.tab | 18 +-
demo/match.R | 101 ++++++++--
man/sim.Rd | 5 +-
24 files changed, 842 insertions(+), 324 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 1693876..3cca3b3 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: Zelig
-Version: 2.4-1
-Date: 2005-08-15
+Version: 2.4-4
+Date: 2005-09-29
Title: Zelig: Everyone's Statistical Software
Author: Kosuke Imai <kimai at Princeton.Edu>,
Gary King <king at harvard.edu>,
@@ -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 Aug 16 11:45:05 2005; olau
+Packaged: Thu Sep 29 23:01:51 2005; olau
diff --git a/NAMESPACE b/NAMESPACE
index f5b3847..7b21a15 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -13,7 +13,7 @@ importFrom(stats, quasipoisson)
export( dims,
gsource,
- help.zelig,
+ help.zelig,
plot.ci,
plot.zelig,
rocplot,
diff --git a/R/MULTIPLE_NEW/make.parameters.R b/R/MULTIPLE_NEW/make.parameters.R
new file mode 100644
index 0000000..f44349e
--- /dev/null
+++ b/R/MULTIPLE_NEW/make.parameters.R
@@ -0,0 +1,39 @@
+make.parameters <- function(terms, shape = "vector", ancillary = TRUE) {
+ 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")
+ 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
+ vars <- unique(unlist(labs))
+ pars <- matrix(NA, ncol = length(syst), nrow = length(vars))
+ colnames(pars) <- eqns[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 = ":")
+ }
+ if (!is.logical(const)) {
+ 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
+}
diff --git a/R/MULTIPLE_NEW/model.frame.list.R b/R/MULTIPLE_NEW/model.frame.list.R
new file mode 100644
index 0000000..2782c94
--- /dev/null
+++ b/R/MULTIPLE_NEW/model.frame.list.R
@@ -0,0 +1,30 @@
+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
new file mode 100644
index 0000000..ede51c9
--- /dev/null
+++ b/R/MULTIPLE_NEW/model.matrix.list.R
@@ -0,0 +1,171 @@
+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
new file mode 100644
index 0000000..1f21f48
--- /dev/null
+++ b/R/MULTIPLE_NEW/parse.formula.R
@@ -0,0 +1,90 @@
+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
new file mode 100644
index 0000000..3bb1a9c
--- /dev/null
+++ b/R/MULTIPLE_NEW/parse.par.R
@@ -0,0 +1,74 @@
+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/put.start.R b/R/MULTIPLE_NEW/put.start.R
new file mode 100644
index 0000000..d59ac03
--- /dev/null
+++ b/R/MULTIPLE_NEW/put.start.R
@@ -0,0 +1,18 @@
+put.start <- function(start.val, value, terms, eqn) {
+ if (!any(class(terms) == "list"))
+ 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 = ":")
+ idx[idx == var] <- paste(idx[idx == var], eqns, collapse = ":")
+ }
+ }
+ par.id <- NULL
+ for (vars in eqn)
+ par.id <- c(par.id, grep(vars, idx))
+ start.val[par.id] <- value
+ start.val
+}
diff --git a/R/MULTIPLE_NEW/set.start.R b/R/MULTIPLE_NEW/set.start.R
new file mode 100644
index 0000000..18b370b
--- /dev/null
+++ b/R/MULTIPLE_NEW/set.start.R
@@ -0,0 +1,15 @@
+set.start <- function(start.val, terms) {
+ if (any(class(terms) == "list"))
+ labs <- make.parameters(terms = terms, shape = "vector", ancillary = TRUE)
+ else
+ labs <- attr(terms, "term.labels")
+ if (is.null(start.val))
+ start.val <- rep(0, length(labs))
+ else {
+ if (length(start.val) != length(labs))
+ stop(paste("length of 'start.val' does not equal number of model parameters = ",
+ length(labs), ".", sep = ""))
+ }
+ names(start.val) <- labs
+ start.val
+}
diff --git a/R/MULTIPLE_NEW/tag.R b/R/MULTIPLE_NEW/tag.R
new file mode 100644
index 0000000..1f88071
--- /dev/null
+++ b/R/MULTIPLE_NEW/tag.R
@@ -0,0 +1,5 @@
+tag <- function (x,y){
+return(x)
+
+}
+
diff --git a/R/MULTIPLE_NEW/terms.list.R b/R/MULTIPLE_NEW/terms.list.R
new file mode 100644
index 0000000..46abcb7
--- /dev/null
+++ b/R/MULTIPLE_NEW/terms.list.R
@@ -0,0 +1,113 @@
+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
new file mode 100644
index 0000000..fb815d4
--- /dev/null
+++ b/R/MULTIPLE_NEW/toBuildFormula.R
@@ -0,0 +1,17 @@
+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/generate.start.R b/R/generate.start.R
deleted file mode 100644
index b509bf0..0000000
--- a/R/generate.start.R
+++ /dev/null
@@ -1,9 +0,0 @@
-generate.start <- function(start.val, X, ancillary = NULL) {
- if (is.null(start.val))
- start.val <- rep(0, ncol(X) + length(ancillary))
- if (!is.null(ancillary))
- names(start.val) <- c(colnames(X), ancillary)
- else
- names(start.val) <- c(colnames(X))
- start.val
-}
diff --git a/R/parse.par.R b/R/parse.par.R
deleted file mode 100644
index 15dea51..0000000
--- a/R/parse.par.R
+++ /dev/null
@@ -1,7 +0,0 @@
-parse.par <- function(par, idx, name) {
- if (name %in% idx)
- return(par[which(idx == name)])
- else
- return(par[which(idx != name)])
-}
-
diff --git a/R/print.summary.MCMCZelig.R b/R/print.summary.MCMCZelig.R
index 5c0e4c0..ebd035c 100644
--- a/R/print.summary.MCMCZelig.R
+++ b/R/print.summary.MCMCZelig.R
@@ -8,6 +8,6 @@ print.summary.MCMCZelig <- function(x, digits=max(3, getOption("digits") -
cat("Sample size per chain =", (x$end -
x$start)/x$thin + 1, "\n")
cat("\n", "Mean, standard deviation, and quantiles for marginal posterior distributions.", "\n")
- print.matrix(round(x$summary, digits=digits))
+ print(round(x$summary, digits=digits))
cat("\n")
}
diff --git a/R/print.summary.zelig.R b/R/print.summary.zelig.R
index 629cde8..815e1d4 100644
--- a/R/print.summary.zelig.R
+++ b/R/print.summary.zelig.R
@@ -42,7 +42,7 @@ print.summary.zelig <- function(x, digits=getOption("digits"),
rownames(tmp[,,j]) <- 1:nrow(tmp[,,j])
if (!is.null(names(tmp[,,j])))
names(tmp[,,j]) <- NULL
- print.matrix(tmp[,,j], digits=digits, ...)
+ print(tmp[,,j], digits=digits, ...)
}
}
else {
diff --git a/R/qi.MCMCZelig.R b/R/qi.MCMCZelig.R
index 42cedbd..c429046 100644
--- a/R/qi.MCMCZelig.R
+++ b/R/qi.MCMCZelig.R
@@ -1,16 +1,11 @@
-
qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
-
model <- object$zelig
qi <- list()
check <- FALSE
-
if (model %in% c("logit.bayes", "probit.bayes", "oprobit.bayes", "mlogit.bayes"))
check <- TRUE
-
if (model %in% c("logit.bayes","probit.bayes", "normal.bayes",
"poisson.bayes","tobit.bayes")) {
-
if (model == "logit.bayes") {
coef <- object$coefficients
eta <- coef %*% t(x)
@@ -18,12 +13,10 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
dimnames(pr) <- dimnames(ev) <- dimnames(eta)
ev <- 1/(1+exp(-eta))
for (i in 1:ncol(ev))
- pr[,i] <- as.character(rbinom(length(ev[,i]), 1, ev[,i]))
-
+ pr[,i] <- as.character(rbinom(length(ev[,i]), 1, ev[,i]))
qi$ev <- ev
qi$pr <- pr
- qi.name <- list(ev = "Expected Values: E(Y|X)", pr="Predicted
- Values: Y|X")
+ qi.name <- list(ev = "Expected Values: E(Y|X)", pr = "Predicted Values: Y|X")
}
else if (model == "probit.bayes") {
coef <- object$coefficients
@@ -35,17 +28,16 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
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")
+ qi.name <- list(ev = "Expected Values: E(Y|X)", pr = "Predicted Values: Y|X")
}
else if (model =="normal.bayes") {
- coef <- object$coefficients[,1:(ncol(object$coefficients)-1)]
- eta <- coef %*% t(x)
- ev <- matrix(NA, nrow = nrow(eta), ncol = ncol(eta))
- dimnames(ev) <- dimnames(eta)
- ev <- eta
- qi$ev <- ev
- qi.name <- list(ev = "Expected Values: E(Y|X)")
+ nvar <- ncol(object$coefficients)
+ coef <- object$coefficients[,1:(nvar-1)]
+ qi$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")
+
}
else if (model =="tobit.bayes") {
coef <- object$coefficients[,1:(ncol(object$coefficients)-1)]
@@ -54,20 +46,15 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
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))
+ ##cev <- eta + sig*(dnorm(L1)-dnorm(L2))/(pnorm(L2)-pnorm(L1))
if (object$below==-Inf) temp1<-0
- else temp1 <- pnorm(L1)*object$below
+ 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 temp2 <- (1-pnorm(L2))*object$above
+ qi$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") {
coef <- object$coefficients
@@ -79,92 +66,64 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
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")
-
+ qi.name <- list(ev = "Expected Values: E(Y|X)",
+ pr="Predicted Values: Y|X")
}
-
-
if (!is.null(x1)) {
- eta1 <- coef %*% t(x1)
+ eta1 <- coef %*% t(x1)
if (model == "logit.bayes") {
ev1 <- 1/(1+exp(-eta1))
- rr <-ev1/ev
- fd <-ev1-ev
-
- qi$fd <- fd
- qi$rr <- rr
-
+ qi$fd <- ev1 - ev
+ qi$rr <- ev1 / ev
qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
-
qi.name$rr <- "Risk Ratios: P(Y=1|X1)/P(Y=1|X)"
}
else if (model == "probit.bayes") {
ev1 <- pnorm(eta1)
- rr <-ev1/ev
- fd <-ev1-ev
-
- qi$fd <- fd
- qi$rr <- rr
-
+ qi$rr <-ev1/ev
+ qi$fd <-ev1-ev
qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
-
qi.name$rr <- "Risk Ratios: P(Y=1|X1)/P(Y=1|X)"
}
else if (model == "normal.bayes") {
ev1 <- eta1
- fd <-ev1-ev
-
- qi$fd <- fd
-
+ qi$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))
-
+ ##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$fd <-ev1-ev
qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
}
else if (model == "poisson.bayes") {
ev1 <- exp(eta1)
- fd <-ev1-ev
-
- qi$fd <- fd
-
+ qi$fd <- exp(eta1) - ev
qi.name$fd <- "First Differences in Expected Values: E(Y|X1)-E(Y|X)"
-
}
- }
-
+ }
if (!is.null(y)) {
- yvar <- matrix(rep(y, nrow(simpar)), nrow = nrow(simpar), byrow = TRUE)
- tmp.ev <- yvar - qi$ev
- if (check)
- tmp.pr <- yvar - as.integer(qi$pr)
- else
- tmp.pr <- yvar - qi$pr
- qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(simpar))
- qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
- if (model %in% c("logit", "probit", "poisson"))
- {
+ yvar <- matrix(rep(y, nrow(simpar)), nrow = nrow(simpar), byrow
+ = TRUE)
+ tmp.ev <- yvar - qi$ev
+ if (check)
+ tmp.pr <- yvar - as.integer(qi$pr)
+ else
+ tmp.pr <- yvar - qi$pr
+ qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(simpar))
+ qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
+ if (model %in% c("logit", "probit", "poisson")) {
qi$ate.pr <- matrix(apply(tmp.pr, 1, mean), nrow = nrow(simpar))
qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
}
- }
-
- list(qi=qi, qi.name=qi.name)
+ }
+ out <- list(qi=qi, qi.name=qi.name)
}
else if ((model =="oprobit.bayes") || (model == "mlogit.bayes")) {
if (model == "oprobit.bayes") {
@@ -172,235 +131,156 @@ qi.MCMCZelig <- function(object, simpar=NULL, x, x1 = NULL, y = NULL, ...) {
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)]
-
-
+ 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]
-
+ ## dimnames(pr)[1] <- dimnames(ev)[1] <- dimnames(eta)[1]
+ ## dimnames(pr)[2] <- dimnames(ev)[3] <- dimnames(eta)[2]
for (j in 1:level)
ev[,j,] <- pnorm(gamma[,j+1]-eta) - pnorm(gamma[,j]-eta)
-
colnames(ev) <- levels(model.response(model.frame(object)))
-
for (j in 1:nrow(pr)) {
- mu <- eta[j,]
-# pr[j,]<-as.character(cut(mu, gamma[j,],
-# labels=as.factor(1:level)))
- pr[j,]<-as.character(cut(mu, gamma[j,], labels=colnames(ev)))
- }
-
-
+ mu <- eta[j,]
+ ## pr[j,]<-as.character(cut(mu, gamma[j,],
+ ## labels=as.factor(1:level)))
+ pr[j,]<-as.character(cut(mu, gamma[j,], labels=colnames(ev)))
+ }
colnames(ev) <- levels(model.response(model.frame(object)))
qi$ev <- ev
qi$pr <- pr
- qi.name <- list(ev = "Expected Values: P(Y=j|X)", pr="Predicted
- Values: Y|X")
+ qi.name <- list(ev = "Expected Values: P(Y=j|X)",
+ pr="Predicted Values: Y|X")
}
else if (model == "mlogit.bayes") {
library(stats)
resp <- model.response(model.frame(object))
level <- length(table(resp))
-
p <- dim(model.matrix(eval(object),data=eval(object$data)))[2]
-
coef <- object$coefficients
-
-
eta <- array(NA, c(nrow(coef),level, nrow(x)))
-
eta[,1,]<-matrix(0, dim(eta)[1],dim(eta)[3])
for (j in 2:level) {
ind <- (1:p)*(level-1)-(level-j)
eta[,j,]<- coef[,ind]%*%t(x)
}
-
eta<-exp(eta)
-
ev <- array(NA, c(nrow(coef), level, nrow(x)))
pr <- matrix(NA, nrow(coef), nrow(x))
-
-
colnames(ev) <- rep(NA, level)
-
- for (k in 1:nrow(x))
- for (j in 1:level)
- {
- ev[,j,k] <- eta[,j,k]/rowSums(eta[,,k])
- }
-
- for (j in 1:level)
- {
- colnames(ev)[j] <- paste("P(Y=", j, ")", sep="")
- }
-
+ for (k in 1:nrow(x)) {
+ for (j in 1:level)
+ ev[,j,k] <- eta[,j,k]/rowSums(eta[,,k])
+ }
+ for (j in 1:level) {
+ colnames(ev)[j] <- paste("P(Y=", j, ")", sep="")
+ }
for (k in 1:nrow(x)) {
probs <- as.matrix(ev[,,k])
temp <- apply(probs, 1, FUN=rmultinom, n=1, size=1)
temp <- as.matrix(t(temp)%*%(1:nrow(temp)))
pr <- apply(temp,2,as.character)
- }
-
+ }
qi$ev <- ev
qi$pr <- pr
- qi.name <- list(ev = "Expected Values: P(Y=j|X)", pr="Predicted
- Values: Y|X")
+ qi.name <- list(ev = "Expected Values: P(Y=j|X)",
+ pr = "Predicted Values: Y|X")
}
-
-
if (!is.null(x1)) {
-
if (model == "oprobit.bayes") {
-
eta1 <- coef %*% t(x1)
-
- ev1 <- array(NA, c(nrow(eta), level, ncol(eta)))
-
- for (j in 1:level)
- ev1[,j,] <- pnorm(gamma[,j+1]-eta1) - pnorm(gamma[,j]-eta1)
-
- rr <-ev1/ev
- fd <-ev1-ev
-
- qi$fd <- fd
- qi$rr <- rr
-
-
+ 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)
+ qi$rr <-ev1/ev
+ qi$fd <-ev1-ev
qi.name$fd <- "First Differences in Expected Values: P(Y=j|X1)-P(Y=j|X)"
-
qi.name$rr <- "Risk Ratios: P(Y=j|X1)/P(Y=j|X)"
-
}
else if (model == "mlogit.bayes") {
-
- eta1 <- array(NA, c(nrow(coef),level, nrow(x1)))
-
- eta1[,1,]<-matrix(0, dim(eta1)[1],dim(eta1)[3])
- for (j in 2:level) {
- ind <- (1:p)*(level-1)-(level-j)
- eta1[,j,]<- coef[,ind]%*%t(x1)
- }
-
- eta1<-exp(eta1)
-
- ev1 <- array(NA, c(nrow(eta1), level, nrow(x1)))
-
-
- for (k in 1:nrow(x))
- for (j in 1:level)
- {
- ev1[,j,k] <- eta1[,j,k]/rowSums(eta1[,,k])
-
- }
-
- rr <-ev1/ev
- fd <-ev1-ev
-
- qi$fd <- fd
- qi$rr <- rr
-
-
+ 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])
+ }
+ qi$rr <-ev1/ev
+ qi$fd <-ev1-ev
qi.name$fd <- "First Differences in Expected Values: P(Y=j|X1)-P(Y=j|X)"
-
qi.name$rr <- "Risk Ratios: P(Y=j|X1)/P(Y=j|X)"
-
- }
+ }
}
-
if (!is.null(y)) {
- yvar <- matrix(rep(y, nrow(simpar)), nrow = nrow(simpar), byrow =
- TRUE)
-
- levels.names<-levels(as.factor(y))
- yvar1<- pr1 <-array(NA, c(nrow(yvar), level, ncol(yvar)))
-
-
- for (j in 1:nrow(yvar))
- {
+ 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))
+ pr1[j,,]<-t(class.ind(as.integer(qi$pr[j,]), levels.names))
else
pr1[j,,]<-t(class.ind(qi$pr[j,],levels.names))
}
-
- tmp.ev <- yvar1 - qi$ev
- tmp.pr <- yvar1 - pr1
- qi$ate.ev <- matrix(apply(tmp.ev, 2, rowMeans), nrow = nrow(simpar))
- qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
- if (model %in% c("oprobit.bayes", "mlogit.bayes"))
- {
+ tmp.ev <- yvar1 - qi$ev
+ tmp.pr <- yvar1 - pr1
+ qi$ate.ev <- matrix(apply(tmp.ev, 2, rowMeans), nrow = nrow(simpar))
+ qi.name$ate.ev <- "Average Treatment Effect: Y - EV"
+ if (model %in% c("oprobit.bayes", "mlogit.bayes", "normal.bayes")) {
qi$ate.pr <- matrix(apply(tmp.pr, 2, rowMeans), nrow = nrow(simpar))
qi.name$ate.pr <- "Average Treatment Effect: Y - PR"
}
- }
-
- list(qi=qi, qi.name=qi.name)
-
-
-
+ }
+ out <- list(qi=qi, qi.name=qi.name)
}
else if (model == "ei.hier" || model == "ei.dynamic") {
if (!any(class(x)=="cond")) stop("set 'cond=TRUE' in setx.\n")
- else
- {
- coef <- object$coefficients
- n <- nrow(x)
- if (is.null(object$N))
- N<-rep(1,n)
- else N <- eval(object$N)
- ev <- array(NA, c(nrow = nrow(coef), 2,2, n))
- pr <- array(NA, c(nrow = nrow(coef), 2,2, n))
- nlen<-length(coef[,1])
- for (j in 1:2) {
- ev[,j,1,] <- t(apply(coef[,((1:n)+(j-1)*n)],
- 1,"*", x[,j])*N)
- ev[,j,2,] <- t(apply((1-coef[,((1:n)+(j-1)*n)]), 1,"*",
- x[,j])*N)
- for (i in 1:n)
- {
- size<-round(x[i,j]*N[i])
- pr[,j,1,i] <-rbinom(prob=coef[,(i+(j-1)*n)], n=nlen, size=size)
-
- pr[,j,2,i] <- x[i,j]*N[i]-pr[,j,1,i]
- }
-
+ else {
+ coef <- object$coefficients
+ n <- nrow(x)
+ if (is.null(object$N))
+ N<-rep(1,n)
+ else N <- eval(object$N)
+ ev <- array(NA, c(nrow = nrow(coef), 2,2, n))
+ pr <- array(NA, c(nrow = nrow(coef), 2,2, n))
+ nlen<-length(coef[,1])
+ for (j in 1:2) {
+ ev[,j,1,] <- t(apply(coef[,((1:n)+(j-1)*n)], 1,"*", x[,j])*N)
+ ev[,j,2,] <- t(apply((1-coef[,((1:n)+(j-1)*n)]), 1,"*", x[,j])*N)
+ for (i in 1:n) {
+ size<-round(x[i,j]*N[i])
+ pr[,j,1,i] <-rbinom(prob=coef[,(i+(j-1)*n)], n=nlen, size=size)
+ pr[,j,2,i] <- x[i,j]*N[i]-pr[,j,1,i]
}
-# dimnames(ev)[[1]] <- dimnames(pr)[[4]] <- 1:nrow(coef)
- dimnames(ev)[[4]] <- dimnames(pr)[[4]] <- rownames(x)
- dimnames(ev)[[2]] <- dimnames(pr)[[2]] <- colnames(x)
- dimnames(ev)[[3]] <- dimnames(pr)[[3]] <- colnames(model.response(object$model))
- class(ev) <- class(pr) <- c("ei", "array")
- qi$ev <- ev
- qi$pr <- pr
- qi.name <- list(ev = "Expected In sample predictions at aggregate
- level", pr = "In sample predictions at aggregate level")
-
}
-
- list(qi=qi, qi.name=qi.name)
- }
- else if ( model %in% c("factor.bayes", "factor.ord",
- "factor.mix", "irt1d", "irtkd"))
- {
- stop("sim procedure not applicable since no explanatory
- variables are involved.\n")
- list(qi=qi)
+ ## dimnames(ev)[[1]] <- dimnames(pr)[[4]] <- 1:nrow(coef)
+ dimnames(ev)[[4]] <- dimnames(pr)[[4]] <- rownames(x)
+ dimnames(ev)[[2]] <- dimnames(pr)[[2]] <- colnames(x)
+ dimnames(ev)[[3]] <- dimnames(pr)[[3]] <- colnames(model.response(object$model))
+ class(ev) <- class(pr) <- c("ei", "array")
+ qi$ev <- ev
+ qi$pr <- pr
+ qi.name <- list(ev = "Expected In sample predictions at aggregate level",
+ pr = "In sample predictions at aggregate level")
}
-
-
-
+ out <- list(qi=qi, qi.name=qi.name)
+ }
+ else if ( model %in% c("factor.bayes", "factor.ord", "factor.mix", "irt1d", "irtkd")) {
+ stop("sim procedure not applicable since no explanatory variables are involved.\n")
+ out <- list(qi=qi)
+ }
+ out
}
diff --git a/R/qi.survreg.R b/R/qi.survreg.R
index f9c78c2..f22259d 100644
--- a/R/qi.survreg.R
+++ b/R/qi.survreg.R
@@ -56,31 +56,36 @@ qi.survreg <- function(object, simpar, x, x1 = NULL, y = NULL) {
qi.name$fd <- "First Differences: E(Y|X1)-E(Y|X)"
}
if (!is.null(y)) {
- tmp <- list(ev = ev$ev[, which(status == 0)], theta = ev$theta[, which(status == 0)])
- y.obs <- matrix(y[status == 1], nrow = nrow(qi$ev),
- ncol = length(y[status == 1]), byrow = TRUE)
- y.imp <- matrix(NA, nrow = nrow(qi$ev), ncol = length(y[status == 0]))
- tmp.scale <- c(matrix(sim.scale, nrow = length(sim.scale),
- ncol = length(y[status == 0])))
- y.imp <- matrix(pr.surv(model, tmp$theta, tmp.scale, tmp$ev),
- nrow = nrow(qi$ev), ncol = length(y[status == 0]))
- y.c <- y[status == 0]
- idx <- t(apply(y.imp, 1, '>=', y.c))
- count <- 1
- while ((sum(idx) < length(idx)) & count < 1001) {
- count <- count + 1
- tmp.idx <- which(!idx, arr.ind = TRUE)
- y.imp[tmp.idx] <- pr.surv(model, tmp$theta[tmp.idx],
- sim.scale[tmp.idx[,1]], tmp$ev[tmp.idx])
- idx[tmp.idx] <- y.imp[tmp.idx] >= y.c[tmp.idx[,2]]
+ if (any(status == 0)) {
+ tmp <- list(ev = ev$ev[, which(status == 0)],
+ theta = ev$theta[, which(status == 0)])
+ y.obs <- matrix(y[status == 1], nrow = nrow(qi$ev),
+ ncol = length(y[status == 1]), byrow = TRUE)
+ y.imp <- matrix(NA, nrow = nrow(qi$ev), ncol = length(y[status == 0]))
+ tmp.scale <- c(matrix(sim.scale, nrow = length(sim.scale),
+ ncol = length(y[status == 0])))
+ y.imp <- matrix(pr.surv(model, tmp$theta, tmp.scale, tmp$ev),
+ nrow = nrow(qi$ev), ncol = length(y[status == 0]))
+ y.c <- y[status == 0]
+ idx <- t(apply(y.imp, 1, '>=', y.c))
+ count <- 1
+ while ((sum(idx) < length(idx)) & count < 1001) {
+ count <- count + 1
+ tmp.idx <- which(!idx, arr.ind = TRUE)
+ y.imp[tmp.idx] <- pr.surv(model, tmp$theta[tmp.idx],
+ sim.scale[tmp.idx[,1]], tmp$ev[tmp.idx])
+ idx[tmp.idx] <- y.imp[tmp.idx] >= y.c[tmp.idx[,2]]
+ }
+ if (count == 1001) {
+ warning(" Maximum number of imputed values (1000) reached for censored Y. \n Using censoring point as observed value, since Pr(Y > Yc | sims) <= 0.001.")
+ y.imp[which(idx == 0, arr.ind = TRUE)] <- y.c[which(idx == 0, arr.ind == TRUE)[,2]]
+ }
+ yvar <- matrix(NA, ncol = length(y), nrow = nrow(qi$ev))
+ yvar[, which(status == 1)] <- y.obs
+ yvar[, which(status == 0)] <- y.imp
}
- if (count == 1001) {
- warning(" Maximum number of imputed values (1000) reached for censored Y. \n Using censoring point as observed value, since Pr(Y > Yc | sims) <= 0.001.")
- y.imp[which(idx == 0, arr.ind = TRUE)] <- y.c[which(idx == 0, arr.ind == TRUE)[,2]]
- }
- yvar <- matrix(NA, ncol = length(y), nrow = nrow(qi$ev))
- yvar[, which(status == 1)] <- y.obs
- yvar[, which(status == 0)] <- y.imp
+ else
+ yvar <- matrix(y, ncol = length(y), nrow = nrow(qi$ev), byrow = TRUE)
tmp.ev <- yvar - qi$ev
tmp.pr <- yvar - qi$pr
qi$ate.ev <- matrix(apply(tmp.ev, 1, mean), nrow = nrow(simpar))
diff --git a/R/setx.default.R b/R/setx.default.R
index ab09ddb..a0ddbda 100644
--- a/R/setx.default.R
+++ b/R/setx.default.R
@@ -58,7 +58,7 @@ setx.default <- function(object, fn = list(numeric = mean, ordered =
odta <- eval(object$call$data, sys.parent())
else
odta <- data
- data <- dta <- model.frame(delete.response(tt), odta)
+ data <- dta <- odta[attributes(model.frame(tt, odta))$row.names,]
vars <- names(dta)
if (!is.null(counter)) {
if (!any(counter == vars))
diff --git a/README b/README
index 970c5db..1880366 100644
--- a/README
+++ b/README
@@ -1,3 +1,11 @@
+2.4-4 (September 29, 2005): Stable release for R 2.0.0-2.2.0. Fixed
+ links for help.zelig().
+
+2.4-3 (September 29, 2005): Stable release for R 2.0.0-2.2.0.
+
+2.4-2 (August 30, 2005): Stable release for R 2.0.0-2.1.1. Fixed bug in
+ setx() related to as.factor() and I(). Streamlined qi.survreg().
+
2.4-1 (August 15, 2005): Stable release for R 2.0.0-2.1.1. Added the
following Bayesian models: factor analysis, mixed factor analysis,
ordinal factor analysis, unidimensional item response theory,
diff --git a/data/MatchIt.url.tab b/data/MatchIt.url.tab
index 9cf2ee9..f647e0c 100644
--- a/data/MatchIt.url.tab
+++ b/data/MatchIt.url.tab
@@ -1,7 +1,4 @@
-matchit http://gking.harvard.edu/matchit/docs/Usage.html
-distance http://gking.harvard.edu/matchit/docs/_TT_distance_TT.html
-matchdef http://gking.harvard.edu/matchit/docs/_TT_matchdef_TT.html
-subclassify http://gking.harvard.edu/matchit/docs/_TT_subclassify_TT.html
-print.matchit http://gking.harvard.edu/matchit/docs/Print.html
-summary.matchit http://gking.harvard.edu/matchit/docs/Summary.html
-plot.matchit http://gking.harvard.edu/matchit/docs/Plot.html
+MatchIt http://gking.harvard.edu/matchit
+matchit http://gking.harvard.edu/matchit/docs/_TT_matchit_TT__Implem.html
+match.data http://gking.harvard.edu/matchit/docs/_TT_match_data_TT__Ext.html
+
diff --git a/data/Zelig.url.tab b/data/Zelig.url.tab
index b6954bd..2a72c0b 100644
--- a/data/Zelig.url.tab
+++ b/data/Zelig.url.tab
@@ -1,7 +1,7 @@
command http://gking.harvard.edu/zelig/docs/Main_Commands.html
commands http://gking.harvard.edu/zelig/docs/Main_Commands.html
-model http://gking.harvard.edu/zelig/docs/Specific_Models.html
-models http://gking.harvard.edu/zelig/docs/Specific_Models.html
+model http://gking.harvard.edu/zelig/docs/Models_Zelig_Can.html
+models http://gking.harvard.edu/zelig/docs/Models_Zelig_Can.html
zelig http://gking.harvard.edu/zelig/docs/_TT_zelig_TT__Estimati.html
setx http://gking.harvard.edu/zelig/docs/_TT_setx_TT__Setting_E.html
sim http://gking.harvard.edu/zelig/docs/_TT_sim_TT__Simulating.html
@@ -20,17 +20,27 @@ repl http://gking.harvard.edu/zelig/docs/_TT_repl_TT__Replicati.html
blogit http://gking.harvard.edu/zelig/docs/_TT_blogit_TT__Bivaria.html
bprobit http://gking.harvard.edu/zelig/docs/_TT_bprobit_TT__Bivari.html
exp http://gking.harvard.edu/zelig/docs/_TT_exp_TT__Exponentia.html
+factor.bayes http://gking.harvard.edu/zelig/docs/_TT_factor_bayes_TT__B.html
+factor.mix http://gking.harvard.edu/zelig/docs/_TT_factor_mix_TT__Mix.html
+factor.ord http://gking.harvard.edu/zelig/docs/_TT_factor_ord_TT__Ord.html
gamma http://gking.harvard.edu/zelig/docs/_TT_gamma_TT__Gamma_Re.html
+irt1d http://gking.harvard.edu/zelig/docs/_TT_irt1d_TT__One_Dime.html
+irtkd http://gking.harvard.edu/zelig/docs/_TT_irtkd_TT__tex2htm.html
ls http://gking.harvard.edu/zelig/docs/_TT_ls_TT__Least_Squar.html
logit http://gking.harvard.edu/zelig/docs/_TT_logit_TT__Logistic.html
+logit.bayes http://gking.harvard.edu/zelig/docs/_TT_logit_bayes_TT__Ba.html
lognorm http://gking.harvard.edu/zelig/docs/_TT_lognorm_TT__Log_No.html
mlogit http://gking.harvard.edu/zelig/docs/_TT_mlogit_TT__Multino.html
-mloglm http://gking.harvard.edu/zelig/docs/_TT_mloglm_TT__Multino.html
+mlogit.bayes http://gking.harvard.edu/zelig/docs/_TT_mlogit_bayes_TT__B.html
negbin http://gking.harvard.edu/zelig/docs/_TT_negbin_TT__Negativ.html
normal http://gking.harvard.edu/zelig/docs/_TT_normal_TT__Normal.html
+normal.bayes http://gking.harvard.edu/zelig/docs/_TT_normal_bayes_TT__B.html
ologit http://gking.harvard.edu/zelig/docs/_TT_ologit_TT__Ordinal.html
oprobit http://gking.harvard.edu/zelig/docs/_TT_oprobit_TT__Ordina.html
+oprobit.bayes http://gking.harvard.edu/zelig/docs/_TT_oprobit_bayes_TT.html
poisson http://gking.harvard.edu/zelig/docs/_TT_poisson_TT__Poisso.html
+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
-weibull http://gking.harvard.edu/zelig/docs/_TT_weibull_TT__Weibul.html
\ No newline at end of file
+weibull http://gking.harvard.edu/zelig/docs/_TT_weibull_TT__Weibul.html
diff --git a/demo/match.R b/demo/match.R
index f6c94f2..036e044 100644
--- a/demo/match.R
+++ b/demo/match.R
@@ -1,45 +1,110 @@
+###
+### Example 1: calculating the average treatment effect for the treated
+###
+
+## load the Lalonde data
library(MatchIt)
data(lalonde)
user.prompt()
-## an example for propensity score matching
-m.out1 <- matchit(treat ~ age + educ + black + hispan + married +
- nodegree + re74 + re75, data = lalonde)
+
+## propensity score matching
+m.out1 <- matchit(treat ~ age + educ + black + hispan + nodegree + married + re74 + re75,
+ method = "nearest", data = lalonde)
user.prompt()
-z.out1 <- zelig(re78 ~ distance + age + educ + black + hispan + married +
- nodegree + re74 + re75,
- data = match.data(m.out1, "control"), model = "ls")
+## fit the linear model to the control group controlling for propensity score and
+## other covariates
+z.out1 <- zelig(re78 ~ age + educ + black + hispan + nodegree + married + re74 + re75 +
+ distance, data = match.data(m.out1, "control"), model = "ls")
user.prompt()
-x.out1 <- setx(z.out1, fn = NULL, data = match.data(m.out1, "treat"), cond = TRUE)
+## set the covariates to the covariates of matched treated units
+## use conditional prediction by setting cond = TRUE.
+x.out1 <- setx(z.out1, data = match.data(m.out1, "treat"), fn = NULL, cond = TRUE)
user.prompt()
+## simulate quantities of interest
s.out1 <- sim(z.out1, x = x.out1)
user.prompt()
+
+## obtain a summary
summary(s.out1)
user.prompt()
-## an example for subclassification
-m.out2 <- matchit(treat ~ age + educ + black + hispan + married +
- nodegree + re74 + re75, data = lalonde,
- replace = FALSE, subclass = 3)
-user.prompt()
-z.out2 <- zelig(re78 ~ distance + age + educ +
- married + nodegree + re74 + re75, data = match.data(m.out2, "control"),
- model="ls", by="subclass")
+###
+### Example 2: calculating the average treatment effect for the entire sample
+###
+
+## fit the linear model to the treatment group controlling for propensity score and
+## other covariates
+z.out2 <- zelig(re78 ~ age + educ + black + hispan + nodegree + married + re74 + re75 +
+ distance, data = match.data(m.out1, "control"), model = "ls")
user.prompt()
-x.out2 <- setx(z.out2, fn = NULL, data = match.data(m.out2, "treat"), cond = TRUE)
+## conducting the simulation procedure for the control group
+x.out2 <- setx(z.out2, data = match.data(m.out1, "control"), fn = NULL, cond = TRUE)
user.prompt()
s.out2 <- sim(z.out2, x = x.out2)
user.prompt()
-summary(s.out2) # overall results
+## Note that Zelig calculates the difference between observed and
+## either predicted or expected values. This means that the treatment
+## effect for the control units is actually the effect of control
+## (observed control outcome minus the imputed outcome under treatment
+## from the model). Hence, to combine treatment effects just reverse
+## the signs of the estimated treatment effect of controls.
+ate.all <- c(s.out1$qi$ate.ev, -s.out2$qi$ate.ev)
user.prompt()
-summary(s.out2, subset = 2) # subclass 2
+## some summaries
+## point estimate
+mean(ate.all)
+user.prompt()
+## standard error
+sd(ate.all)
+user.prompt()
+## 95% confidence interval
+quantile(ate.all, c(0.025, 0.975))
+user.prompt()
+
+
+###
+### Example 3: subclassification
+###
+## subclassification with 4 subclasses
+m.out2 <- matchit(treat ~ age + educ + black + hispan + nodegree + married + re74 + re75,
+ data = lalonde, method = "subclass", subclass = 4)
+user.prompt()
+
+## controlling only for the estimated prpensity score and lagged Y within each subclass
+## one can potentially control for more
+z.out3 <- zelig(re78 ~ re74 + re75 + distance, data = match.data(m.out2, "control"),
+ model = "ls", by = "subclass")
+user.prompt()
+
+## conducting simulations
+x.out3 <- setx(z.out3, data = match.data(m.out2, "treat"), fn = NULL, cond = TRUE)
+user.prompt()
+
+## for the demonstration purpose, we set the number of simulations to be 100
+s.out3 <- sim(z.out3, x = x.out3, num = 100)
+user.prompt()
+
+## overall results
+summary(s.out3)
+user.prompt()
+
+## summary for each subclass
+summary(s.out3, subset = 1)
+user.prompt()
+
+summary(s.out3, subset = 2)
+user.prompt()
+
+summary(s.out3, subset = 3)
+user.prompt()
diff --git a/man/sim.Rd b/man/sim.Rd
index 0e1c4c6..7d28814 100644
--- a/man/sim.Rd
+++ b/man/sim.Rd
@@ -17,7 +17,7 @@
\usage{
s.out <- sim(object, x, x1 = NULL, num = c(1000, 100), prev = NULL,
- bootstrap = FALSE, bootfn = NULL, cond.data = NULL, \dots)
+ bootstrap = FALSE, bootfn = NULL, \dots)
}
\arguments{
@@ -48,9 +48,6 @@ s.out <- sim(object, x, x1 = NULL, num = c(1000, 100), prev = NULL,
bootstrap methods include sampling the residuals rather than the
observations, weighted sampling, and parametric bootstrapping.
(Not available for conditional prediction.) }
- \item{cond.data}{a data frame, identical to the \code{data}
- argument in \code{\link{setx}}. Required for conditional prediction
- with the exponential, Weibull, and lognormal models. }
\item{\dots}{additional optional arguments passed to
\code{boot}. }
}
--
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