Nothing
#
# multinomRob
#
# Walter R. Mebane, Jr.
# University of Michigan
# http://www-personal.umich.edu/~wmebane
# <wmebane@umich.edu>
#
# Jasjeet Singh Sekhon
# UC Berkeley
# http://sekhon.polisci.berkeley.edu
# sekhon@berkeley.edu
#
# $Id: multinomRob.R,v 1.16 2005/09/27 08:04:06 wrm1 Exp $
######################################################
## multinomRob: genoud and Gauss-Newton estimation ##
######################################################
multinomRob <-
function(model,
data,
starting.values=NULL,
equality=NULL, # list of lists of parameter equality constraints
genoud.parms = NULL, ## can be NULL, single, or list
print.level = 0,
iter = FALSE, # should we iterate
maxiter = 10, # maximum number of iterations before we stop
multinom.t=1, # 0=no, 1=yes, 2=force should we do multinom-t for starting values
multinom.t.df=NA, #if set, the multivariate-t function is FORCED to use this DF
MLEonly=FALSE
)
{
#check input
if (print.level < 0)
print.level <- 0;
if(MLEonly!=0 & MLEonly!=1)
MLEonly <- FALSE
if (iter!=TRUE & iter!=FALSE)
{
stop("multinomRob(): illegal input for variable iter")
}
if (maxiter < 0)
{
stop("multinomRob(): illegal input for variable maxiter")
}
if (multinom.t!=0 & multinom.t!=1 & multinom.t!=2)
{
stop("multinomRob(): illegal input for variable multinom.t")
}
if (!is.na(multinom.t.df) & multinom.t.df <= 0)
{
stop("multinomRob(): illegal input for variable multinom.t.df")
}
converged.test <- function(old,new,tol=1e-8) {
abs((old-new)/old) < tol;
}
genoud.fun <- function(z) {
fit.multinomial.C.lqd2(z,X,Y,Ypos,xvec,tvec,ncats,nvars,nvars.unique,obs,TotalY)
}
## ----------------- BEGIN: PARSING / MODEL BUILDING AUTOMATION ------------------
## ---- DATA X, Y and Ypos ----
data.all <- get.xy(model, data, print.level=print.level);
Y <- data.all$Y;
ncats <- dim(Y)[2];
obs <- dim(Y)[1];
Ypos <- data.all$ypos; # element is TRUE if y>=0, FALSE if y<0
Y[!Ypos] <- 0; # set count values to be ignored to zero, for convenience later
X <- data.all$X;
nvars <- dim(X)[2];
## ---- LABELS AND CONTROLS/FLAGS ----
xvec <- matrix(0,nrow=nvars,ncol=ncats)
colnames(xvec) <- choice.labels <- data.all$ynames;
xvar.labels <- rep( "", nvars)
for (i in 1:ncats) {
## only look at named columns
for (j in 1:nvars) {
## and keep only non-padding/zeroed names
ssplit <- ifelse( i==1, "", "/")
if (data.all$xlengths[i] >= j) {
xvar.labels[j] <-
paste( xvar.labels[j], ssplit, data.all$xnames[[i]][j], sep="")
xvec[j,i] <- 1
} else {
xvar.labels[j] <- paste( xvar.labels[j], ssplit, "NA", sep="")
}
}
}
rownames(xvec) <- xvar.labels
# parameter equality constraints
if (!is.null(equality)) {
eqlen <- length(equality);
ieq <- 1;
for (i in 1:eqlen) {
ieq <- ieq + 1;
xynames <- get.xynames(equality[[i]], data);
nynames <- length(xynames$ynames);
yidx <- match(xynames$ynames, data.all$ynames) ;
if (any(is.na(yidx))) {
if (print.level >= 0) {
print("response used in equality constraints is not in the model");
cat("responses in model:", data.all$ynames, "\n");
cat("responses in equality:", xynames$ynames, "\n");
print("terminating multinomRob")
}
return(NULL);
}
xidx <- list();
xidxlen <- 0;
for (j in 1:nynames) {
xidx[[j]] <- match(xynames$xnames[[j]], data.all$xnames[[ yidx[j] ]]);
if (any(is.na(xidx[[j]]))) {
if (print.level >= 0) {
print("regressor used in equality constraints is not in the model");
cat("response:", xynames$ynames[j]);
cat("regressors in model:", data.all$xnames[[ yidx[j] ]], "\n");
cat("regressors in equality:", xynames$xnames, "\n");
print("terminating multinomRob")
}
return(NULL);
}
xidxlen <- xidxlen + length(xidx[[j]]);
}
xyvmatrix <- matrix(0, 3, xidxlen);
kk <- 0;
for (j in 1:nynames) {
for (k in 1:length(xidx[[j]])) {
kk <- kk + 1;
xyvmatrix[1,kk] <- xidx[[j]][k] ;
xyvmatrix[2,kk] <- yidx[j] ;
xyvmatrix[3,kk] <- xvec[xidx[[j]][k], yidx[j]] ;
}
}
if (any(xyvmatrix[3,] != 1)) { # need to consolidate equality constraints
if (any(xyvmatrix[3,1] != xyvmatrix[3,])) {
# not a subset of previous constraints
xyvu <- unique(xyvmatrix[3,]);
xyvu <- xyvu[xyvu>1];
xyvmin <- min(xyvu); # code to consolidate on
for (k in 1:length(xyvu)) {
xvec[xvec == xyvu[k]] <- xyvmin;
xyvmatrix <- xyvmatrix[,xyvmatrix[3,] != xyvu[k]];
}
if (dim(xyvmatrix)[2] > 0) { # parameters not previously constrained
for (k in 1:(dim(xyvmatrix)[2])) {
xvec[xyvmatrix[1,k], xyvmatrix[2,k]] <- xyvmin;
}
}
}
}
else { # consolidation not needed
for (k in 1:(dim(xyvmatrix)[2])) {
xvec[xyvmatrix[1,k], xyvmatrix[2,k]] <- ieq;
}
}
}
if (print.level > 0) {
cat("\nEquality constraints among parameters (after consolidation):\n");
nidxvals <- length(idxvals <- sort(unique(xvec[xvec>1])));
for (i in 1:nidxvals) {
cat("Equality constrained set", i, "\n");
for (j in 1:nvars) for (k in 1:ncats) {
if (xvec[j,k]==idxvals[i]) {
cat("outcome", data.all$ynames[k], "regressor",
data.all$xnames[[k]][j], "\n");
}
}
}
}
}
if (print.level > 0) {
cat("\nYour Model (xvec):\n")
print(xvec)
}
TotalY <- apply( Y, 1, sum)
tvec <- xvec;
nvars.unique <- sum(xvec == 1) + length(unique(xvec[xvec>1]));
# nvars.total <- nvars
## ------------------- END: PARSING / MODEL BUILDING AUTOMATION ------------------
#let's create jacstack
jacstack <- jacstack.function(X,nvars.unique,xvec)
# check for regressors with a distinct value at only one observation
jsingle <- jacstack.singles(jacstack);
if (any(jsingle) & print.level > 0) {
cat("\n");
print("multinomRob: WARNING. Limited regressor variation...")
print("WARNING. ... A regressor has a distinct value for only one observation.")
print("WARNING. ... I'm using a modified estimation algorithm (i.e., preventing LQD")
print("WARNING. ... from modifying starting values for the affected parameters).")
xsingle <- tvec != tvec;
xsingle <- mnl.xvec.mapping(forward=FALSE,xvec,xsingle,jsingle,ncats,nvars);
print("WARNING. ... Affected parameters are TRUE in the following table.")
cat("\n");
print(xsingle);
cat("\n");
}
#load up genoudParms
if(is.null(genoud.parms))
genoud.parms <- list(Domains=NULL);
genoud.parms <- genoudParms(genoud.parms)
# initialize with terrible fit values
mnl.fit <- 9999999999;
multinomT.fit <- 9999999999;
starting.fit <- 9999999999;
multinomT.foo <- NULL
if (is.null(starting.values)) {
starting.values <- vector(mode="numeric", length=nvars.unique);
if (print.level > 0)
cat("\n\nmultinomRob(): Grouped MNL Estimation\n");
Yp <- Y/TotalY
Yp[,ncats] <- 1-apply(as.matrix(Yp[,1:(ncats-1)]),1,sum);
mnl1 <- multinomMLE(Y=Y, Ypos=Ypos, Xarray=X, xvec=xvec, jacstack=jacstack,
xvar.labels=xvar.labels, choice.labels=choice.labels,
MLEonly=MLEonly,
print.level=print.level);
if(!MLEonly)
{
starting.values <- mnl1$coeffvec
mnl.fit <-
fit.multinomial.C.lqd2(starting.values,X,Y,Ypos,xvec,tvec,ncats,nvars,
nvars.unique,obs,TotalY);
}
if (print.level > 0) {
if(!MLEonly)
cat("MNL LQD Fit:",mnl.fit,"\n");
cat("MNL Estimates:\n");
print(mnl1$coefficients);
cat("\n");
cat("MNL SEs:\n");
print(mnl1$se);
cat("\n");
}
if(MLEonly)
{
z <- list(coefficients=mnl1$coefficients,
se=mnl1$se,
mnl=mnl1,
fitted.values=mnl1$fitted.prob,
deviance=as.double(mnl1$GNlist$LLvals)*2,
value=as.double(mnl1$GNlist$LLvals),
error=mnl1$error,
type="MLEonly"
)
class(z) <- "multinomRob"
return(z)
}
starting.values.user <- starting.values
tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,starting.values,
ncats,nvars)
starting.fit <-
fit.multinomial.C.lqd2(starting.values,X,Y,Ypos,xvec,tvec,ncats,nvars,
nvars.unique,obs,TotalY)
use.mnl <- TRUE;
if (multinom.t > 0 & all(Ypos)) {
if (print.level > 0) {
cat("\n\nmultinomRob(): Calculating multinomial-t starting values.\n");
}
Sres.raw <- res.std(Y,TotalY,mnl1$fitted.prob);
multinomT.startvalues <- list()
multinomT.startvalues$beta <- starting.values;
multinomT.startvalues$Omega <- var(Sres.raw);
#multinomT.startvalues$Omega <- var(mnl1$res[,1:(ncats-1)])/(obs-nvars.unique);
multinomT.startvalues$df <- 10;
multinomT.foo <- multinomT(Yp=Yp, Xarray=X, xvec=xvec, jacstack=jacstack,
start=multinomT.startvalues, nobsvec=TotalY,
fixed.df=multinom.t.df);
for (ii in 0:6) {
if(is.null(multinomT.foo$se$beta[1]) & is.na(multinom.t.df)) {
if (print.level > 2)
cat("Fixing Bad Estimate:",ii,"\n")
multinomT.foo <- multinomT(Yp=Yp, Xarray=X, xvec=xvec, jacstack=jacstack,
start=multinomT.startvalues, nobsvec=TotalY,
fixed.df=(10^ii-.5));
if (print.level > 2)
cat("DF: ",multinomT.foo$par$df,"\n");
multinomT.startvalues$beta <- multinomT.foo$par$beta
multinomT.startvalues$Omega <- multinomT.foo$par$Omega
multinomT.startvalues$df <- multinomT.foo$par$df
multinomT.foo <- multinomT(Yp=Yp, Xarray=X, xvec=xvec, jacstack=jacstack,
start=multinomT.startvalues, nobsvec=TotalY);
} else if (is.null(multinomT.foo$se$beta[1]) & !is.na(multinom.t.df)) {
if (print.level > 1)
{
cat("WARNING: Multinom-T SEs are null, but we are using multinomT point estimates for starting values anyways because you explicitly sepecifed a DF\n")
}
use.mnl <- FALSE;
starting.values <- multinomT.foo$par$beta;
multinomT.fit <-
fit.multinomial.C.lqd2(starting.values,X,Y,Ypos,xvec,tvec,ncats,nvars,
nvars.unique,obs,TotalY);
if (print.level > 1)
{
cat("Multinom-T LQD Fit (step ",ii,"):",multinomT.fit,"\n")
beta.print <-
mnl.xvec.mapping(forward=FALSE,xvec,tvec,multinomT.foo$par$beta,ncats,nvars)
cat("Multinom-T Beta Estimates (step ",ii,"):\n")
print( beta.print )
cat("Multinom-T SEs are not defined\n")
}
break;
} else {
use.mnl <- FALSE;
starting.values <- multinomT.foo$par$beta;
multinomT.fit <-
fit.multinomial.C.lqd2(starting.values,X,Y,Ypos,xvec,tvec,ncats,nvars,
nvars.unique,obs,TotalY);
if (print.level > 1)
{
cat("Multinom-T LQD Fit (step ",ii,"):",multinomT.fit,"\n")
beta.print <-
mnl.xvec.mapping(forward=FALSE,xvec,tvec,multinomT.foo$par$beta,ncats,nvars)
se.print <-
mnl.xvec.mapping(forward=FALSE,xvec,tvec,multinomT.foo$se$beta,ncats,nvars)
cat("Multinom-T Beta Estimates (step ",ii,"):\n")
print( beta.print )
cat("Multinom-T Beta SEs (step ",ii,"):\n")
print( se.print )
cat("Multinom-T Omega Estimates (step ",ii,"):\n")
print( multinomT.foo$par$Omega )
cat("Multinom-T DF (step ",ii,"):", multinomT.foo$par$df,"\n\n")
}
break;
}
}#ii
}#end of multinom.t
if (use.mnl==FALSE & !is.null(starting.values) & multinom.t < 2) {
if (multinomT.fit >= mnl.fit)
use.mnl <- TRUE;
}
if(use.mnl==TRUE) {
if (print.level > 0)
cat("multinomRob(): Using grouped MNL estimates as starting values.\n")
starting.values <- vector(mode="numeric", length=nvars.unique);
starting.values <-
mnl.xvec.mapping(forward=TRUE,xvec,mnl1$coefficients,starting.values,
ncats,nvars);
} else {
if (print.level > 0) {
cat("multinomRob(): Using multinomial-t estimates as starting values.\n")
cat("Multinom-T LQD Fit:",multinomT.fit,"\n")
beta.print <-
mnl.xvec.mapping(forward=FALSE,xvec,tvec,multinomT.foo$par$beta,ncats,nvars)
cat("Multinom-T Estimates:\n")
print( beta.print )
cat("Multinom-T DF:", multinomT.foo$par$df,"\n\n")
}
starting.values <- multinomT.foo$par$beta
}
}#is.null(starting.values)
tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,starting.values,ncats,nvars);
fit.starting <-
fit.multinomial.C.lqd2(starting.values,X,Y,Ypos,xvec,tvec,ncats,nvars,
nvars.unique,obs,TotalY)
if (print.level > 2) {
cat("multinomRob(): Starting Values \n")
print(tvec)
cat("multinomRob(): starting fit =",fit.starting,"\n");
}
genoud.out <- genoudRob(genoud.fun,
nvars=nvars.unique,
starting.values = starting.values,
as.list(genoud.parms))
s0 <- genoud.out$value;
lqd.beta.vector <- genoud.out$par;
tvec <- xvec;
tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,lqd.beta.vector,ncats,nvars);
if (print.level>0) {
cat("\n");
cat("\nLQD Results:\n");
print(tvec);
cat("\nLQD sigma:",s0,"\n\n")
}
#############################################
## TANH #
#############################################
cat("\n(multinomTanh):\n");
# modify LQD results for single unique value regressors
tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,
ifelse(jsingle, starting.values, lqd.beta.vector),
ncats,nvars);
mout <- multinomTanh(Y=Y, Ypos=Ypos, X=X, jacstack=jacstack,
xvec=xvec, tvec=as.data.frame(tvec),
pop=TotalY,s2=s0^2,
xvar.labels=xvar.labels, choice.labels,
print.level=print.level);
error <- mout$mtanh$error;
if (mout$mtanh$error==0) {
if (print.level > 0) {
cat("Tanh Estimates\n")
print(mout$mtanh$coef)
cat("\nTanh Sandwich SEs\n");
print(mout$mtanh$se);
cat("\n")
cat("TANH sigma:",sqrt(mout$mtanh$tanhsigma2),"\n\n")
} #end print
} else {
if (print.level >= 0) {
cat("WARNING: Tanh returned an error code:",mout$mtanh$error,"\n");
cat("WARNING: Optimization not complete. Starting values are probably not good enough. \n")
if(iter) {
cat("WARNING: Hopefully we will fix this problem during subsequent iterations.\n")
} else {
cat("WARNING: Rerun with new different starting value or with the iteration option left on\n")
}
}
}
# If the fit value gets worse, we stop. Which is *NOT*
# what the reference file does. The reference file
# keeps going (even if the new value is worse until there is
# stability.
s.count <- 0;
if (iter)
{
s0.old <- s0;
s0.tanh.old <- sqrt(mout$mtanh$tanhsigma2);
tolerance <- genoud.parms$solution.tolerance;
#let's shrink our domains by a half for the iteration stuff because we are
# going to assume that the tanh starting values are not that bad.
genoud.parms$scale.domains <- genoud.parms$scale.domains/2;
mout.old <- mout;
}#iter
while(iter)
{
s.count <- s.count+1;
if (print.level > 0) {
cat("\n**********************************\n\n");
cat("Iteration:",s.count,"\n")
}
if (!any(is.na(mout.old$mtanh$coeffvec))) {
LQDfit.tanhcoefs <-
fit.multinomial.C.lqd2(mout.old$mtanh$coeffvec,X,Y,Ypos,xvec,tvec,
ncats,nvars,nvars.unique,obs,TotalY);
} else {
LQDfit.tanhcoefs <- 9999999999;
}
best <- order(c(starting.fit, mnl.fit, multinomT.fit, LQDfit.tanhcoefs))[1]
if (s.count > 1 && !any(is.na(mout.old$mtanh$coeffvec))) {
if (print.level > 0)
cat("using LQD sigma provided by Tanh.\n")
starting.values <- mout.old$mtanh$coeffvec
} else if (best==1) {
if (print.level > 0)
cat("using best non-genoud LQD sigma. Provided by user starting values.\n")
starting.values <- starting.values.user
} else if (best==2) {
if (print.level > 0)
cat("using best non-genoud LQD sigma. Provided by ML multinomial.\n")
starting.values <- mnl1$coeffvec
} else if (best==3) {
if (print.level > 0)
cat("using best non-genoud LQD sigma. Provided by multinomial-t.\n")
starting.values <- multinomT.foo$par$beta;
} else if (best==4) {
if (print.level > 0)
cat("using best non-genoud LQD sigma. Provided by Tanh.\n")
starting.values <- mout.old$mtanh$coeffvec
}
genoud.out <- genoudRob(genoud.fun,
nvars=nvars.unique,
starting.values = starting.values,
as.list(genoud.parms))
s0 <- genoud.out$value;
lqd.beta.vector <- genoud.out$par;
tvec <- xvec;
tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,lqd.beta.vector,
ncats,nvars);
if (print.level>0) {
cat("\n(iter",s.count,"): LQD Results:\n");
print(tvec);
cat("\n(iter",s.count,"): LQD sigma:",s0,"\n\n")
}
if (s0 > (s0.old+tolerance))
{
if (print.level > 0)
{
cat("WARNING: New fit is worse. Stopping.\n");
cat("It may be a good idea to restart with a larger GENOUD population.\n");
} #end of print.level
error <- -1;
iter <- FALSE;
break;
}
if (!is.na(s0.tanh.old) && converged.test(s0.old, s0, tolerance))
{
if (print.level > 0)
{
cat("\n**********************************\n")
cat("CONVERGED\n");
if (print.level > 0) {
cat("(Converged): Tanh Estimates\n")
print(mout$mtanh$coef)
cat("\n(Converged): Tanh Sandwich SEs\n");
print(mout$mtanh$se);
cat("\n")
cat("(Converged): TANH sigma:",sqrt(mout$mtanh$tanhsigma2),"\n\n")
s0.tanh <- sqrt(mout$mtanh$tanhsigma2);
} #end print
} #end of print.level
error <- 0;
iter <- FALSE;
break;
} else {
#new s0 is smaller
# modify LQD results for single unique value regressors
tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,
ifelse(jsingle, starting.values, lqd.beta.vector),
ncats,nvars);
mout <- multinomTanh(Y=Y, Ypos=Ypos, X=X, jacstack=jacstack,
xvec=xvec, tvec=as.data.frame(tvec),
pop=TotalY,s2=s0^2,
xvar.labels=xvar.labels, choice.labels,
print.level=print.level);
error <- mout$mtanh$error;
if (mout$mtanh$error==0)
{
if (print.level > 0) {
cat("(iter",s.count,"): Tanh Estimates\n")
print(mout$mtanh$coef)
cat("\n(iter",s.count,"): Tanh Sandwich SEs\n");
print(mout$mtanh$se);
cat("\n")
cat("(iter",s.count,"): TANH sigma:",sqrt(mout$mtanh$tanhsigma2),"\n\n")
s0.tanh <- sqrt(mout$mtanh$tanhsigma2);
mout.old <- mout;
s0.tanh.old <- s0.tanh;
} #end print
} else {
if (print.level > 0)
cat("(iter",s.count,"): tanh returns an error:",mout$mtanh$error,"\n");
}
} #end of else
if (s.count==maxiter)
{
if (print.level > 0)
cat("\nWARNING: Maximum number of iterations reached. Results have NOT converged\n\n");
#iter failed
error <- -1;
iter <- FALSE;
break;
}
if (s0 < s0.old)
s0.old <- s0;
}#iter check
if (!exists("mnl1"))
mnl1 <- NULL;
if (!exists("multinomT.foo"))
multinomT.foo <- NULL;
#Rotated Residuals. This will result in a relatively easy to interpret
#vector of residuals. But this residual vector will NOT be a
#consistent set of ortho residuals.
residuals.rotate <- matrix(nrow=obs,ncol=ncats)
if (mout$mtanh$error < 32) {
for (ii in 1:ncats)
{
tindx <- 1:ncats
tindx[1] <- ii;
tindx[ii] <- 1;
YTmp <- Y[,c(tindx)]
YposTmp <- Ypos[,c(tindx)]
XTmp <- X[,,c(tindx), drop=FALSE];
jacstackTmp <- jacstack[,,c(tindx), drop=FALSE]
tvec <- as.matrix(mout$mtanh$coefficients[,c(tindx), drop=FALSE])
foo <- permute(Y=YTmp, Ypos=YposTmp, Xarray=XTmp,
jacstack=jacstackTmp, tvec=tvec,
pop=TotalY, sigma=sqrt(mout$mtanh$disp),
weight=mout$mtanh$w)
residuals.rotate[,ii] <- foo$student[,1];
} #end of ii loop
}
z <- list(coefficients=mout$mtanh$coefficients,
se=mout$mtanh$se,
LQDsigma2=mout$mtanh$dispersion,
TANHsigma2=mout$mtanh$tanhsigma2,
weights=mout$weights,
Hdiag=mout$Hdiag,
prob=mout$mtanh$prob,
residuals.rotate=residuals.rotate,
residuals.student=mout$cr$student,
residuals.standard=mout$cr$standard,
mnl=mnl1,
multinomT = multinomT.foo,
genoud = genoud.out,
mtanh = mout$mtanh,
error = error,
iter = s.count,#the iter number at the solution
type="robust");
class(z) <- "multinomRob"
return(z)
}
#mimicking the definition of summary.default
#summary.default <-
# function(object, ..., digits = max(3, getOption("digits") - 3))
summary.multinomRob <- function(object, ..., digits=3, weights=FALSE)
{
if (class(object) != "multinomRob") {
warning("Object not of class 'multinomRob'")
return(NULL)
}
if(object$type!="MLEonly")
{
out.mtanh <- list()
if (!is.null(object$mtanh) && is.list(object$mtanh) && object$mtanh$error == 0) {
ncats <- NCOL(object$mtanh$coefficients);
nx <- NROW(object$mtanh$coefficients);
# find length of longest regressor label
labmax <- 0;
for (i in 1:ncats) {
for (j in 1:nx) {
labn <- nchar(strsplit(row.names(object$mtanh$se)[j], "/")[[1]][i]) ;
if (labmax < labn) labmax <- labn;
}
}
# generate a blank variable with as many spaces as the longest label
spc <- " ";
blank <- "";
for (i in 1:labmax) blank <- paste(spc, blank, sep="");
for (i in 1:ncats) {
tmp <- as.data.frame(
list("Est" = object$mtanh$coefficients[,i],
"SE Sand" = object$mtanh$se[,i],
# "SE OPG" = object$mtanh$se.opg[,i],
# "SE Hess" = object$mtanh$se.hes[,i],
"t-val Sand" = object$mtanh$coefficients[,i]/object$mtanh$se[,i]))
tmp <- as.matrix(tmp);
rn <- rep(blank, nx);
for (j in 1:nx) {
rnj <- strsplit(row.names(object$mtanh$se)[j], "/")[[1]][i] ;
substr(rn[j], 1, nchar(rnj)) <- rnj;
}
dimnames(tmp)[[1]] <- rn;
out.mtanh[[i]] <- tmp
choice.lables <- labels(object$coefficients)[[2]]
cat("\nChoice",i,":",choice.lables[i],"Estimates and SE:\n")
print(signif(tmp,digits=digits))
cat("\n")
}
cat("\n");
cat("LQD sigma:",sqrt(object$LQDsigma2),"\n")
cat("TANH sigma:",sqrt(object$TANHsigma2),"\n\n")
indx <- object$weights[,2:ncats]==0;
indx[is.na(indx)] <- FALSE;
if (ncats==2) indx <- matrix(indx, length(indx), 1);
w0.obs <- sum(apply(indx,1,any));
cat("Number of Observations:",NROW(object$prob),"\n")
cat("Number of observations with at least one zero weight:",w0.obs,"\n")
w0 <- sum(indx)
cat("Number of zero weights:",w0,"\n")
if (weights) {
cat("TANH: weights\n");
print(object$weights);
}
}
else {
stop("error encountered in tanh result object");
}
} else { #END OF !MLEonly
if (!is.null(object$mnl) && is.list(object$mnl) && object$mnl$error == 0) {
ncats <- NCOL(object$mnl$coefficients);
nx <- NROW(object$mnl$coefficients);
# find length of longest regressor label
labmax <- 0;
for (i in 1:ncats) {
for (j in 1:nx) {
labn <- nchar(strsplit(row.names(object$mnl$se)[j], "/")[[1]][i]) ;
if (labmax < labn) labmax <- labn;
}
}
# generate a blank variable with as many spaces as the longest label
spc <- " ";
blank <- "";
for (i in 1:labmax) blank <- paste(spc, blank, sep="");
for (i in 1:ncats) {
tmp <- as.data.frame(
list("Est" = object$mnl$coefficients[,i],
"SE Sand" = object$mnl$se[,i],
# "SE OPG" = object$mnl$se.opg[,i],
# "SE Hess" = object$mnl$se.hes[,i],
"t-val Sand" = object$mnl$coefficients[,i]/object$mnl$se[,i]))
tmp <- as.matrix(tmp);
rn <- rep(blank, nx);
for (j in 1:nx) {
rnj <- strsplit(row.names(object$mnl$se)[j], "/")[[1]][i] ;
substr(rn[j], 1, nchar(rnj)) <- rnj;
}
dimnames(tmp)[[1]] <- rn;
choice.lables <- labels(object$coefficients)[[2]]
cat("\nChoice",i,":",choice.lables[i],"Estimates and SE:\n")
print(signif(tmp,digits=digits))
cat("\n")
}
cat("\n");
cat("Residual Deviance:",as.double(object$deviance),"\n\n")
} else {
stop("error encountered in MNL result object");
}
}#end of MLEonly
} #end of summary.multinomRob()
#mimicking the definition of plot
plot.multinomRob <- function(x, ...)
{
#this plots the residuals
if (class(x) != "multinomRob") {
warning("Object not of class 'multinomRob'")
return(NULL)
}
# out.mtanh <- list()
if (!is.null(x$mtanh) && is.list(x$mtanh) && x$mtanh$error == 0)
{
ncats <- NCOL(x$mtanh$coefficients);
choice.lables <- labels(x$coefficients)[[2]]
ask.start <- par("ask",no.readonly=TRUE)
for (i in 1:(ncats-1))
{
plot(x$residuals.student[,i],xlab="Observation", ylab="Residual", ...)
title(main=paste(choice.lables[i],"\nStudentized Residuals"));
if (i==1)
par(ask=TRUE);
} #end for
par(ask=ask.start)
}
else {
print("error encountered in tanh result object");
}
#end if
}#end of plot.multinomRob
#Function to create ortho-residuals (with base correction) for the
#other choices. This will result in a relatively easy to interpret
#vector of residuals. But this residual vector will NOT be a
#consistent set of ortho residuals.
#
#Modeled on
#lapo:~/xchg/election/R/multinomial/FL2/base4.permute4.origtanh2.R
#permute.newtanh4.R
#
permute <- function(Y, Ypos, Xarray, jacstack, tvec, pop, sigma, weight)
{
#tvec: tanh estimated values
#from multinomTanh
nobs <- dim(Y)[1]
ncats <- dim(Y)[2]
Hdiag <- robustified.leverage(tvec, Y, Ypos, Xarray, pop, ifelse(weight >0,1,0),jacstack);
# w.Hdiag <- as.data.frame(matrix(c(as.vector(1:nobs),
# signif(weight), signif(Hdiag)),ncol=(ncats-1)+(ncats-1)+1));
# names(w.Hdiag) <-
# c("name",paste("weights:",choice.labels[1:ncats-1],sep=""),
# paste("Hdiag:",choice.labels[1:ncats-1],sep=""));
#cat("mtanh: weights, Hdiag (by choices)\n");
cr <- fn.region.results(tvec, Y, Ypos, Xarray, pop, sigma, Hdiag);
return(list(pred=cr$pred, student=cr$student, standard=cr$standard, Hdiag=Hdiag))
} #end of permute
.onAttach <- function( ... )
{
MatchLib <- dirname(system.file(package = "multinomRob"))
version <- packageDescription("multinomRob", lib.loc = MatchLib)$Version
BuildDate <- packageDescription("multinomRob", lib.loc = MatchLib)$Date
foo <- paste("## \n## multinomRob (Version ", version, ", Build Date: ", BuildDate, ")\n",
"## See http://sekhon.berkeley.edu/robust for additional documentation.\n",
"## Please cite as: Walter R. Mebane, Jr. and Jasjeet S. Sekhon. \"Robust Estimation\n",
"## and Outlier Detection for Overdispersed Multinomial Models of Count Data\".\n",
"## American Journal of Political Science, 48 (April): 391-410. 2004.\n##\n",
sep="")
packageStartupMessage(foo)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.