# R/multinomRob.R In multinomRob: Robust Estimation of Overdispersed Multinomial Regression Models

#### Documented in multinomRob

```#
#  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) {
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");
}

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)
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]]

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)
} #end for
}
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)
}
```

## Try the multinomRob package in your browser

Any scripts or data that you put into this service are public.

multinomRob documentation built on May 2, 2019, 6:54 a.m.