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

#### Defines functions jacstack.functionjacstack.singles

#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon
#  UC Berkeley
#  http://sekhon.polisci.berkeley.edu
#  sekhon@berkeley.edu
#
#  \$Id: datamanip.R,v 1.7 2005/09/27 08:04:06 wrm1 Exp \$
#

#Mapping from xvec to the beta.vector and back again
#forward==TRUE from xvec TO beta.vector
#forward==FALSE from beta.vector TO xvec
mnl.xvec.mapping <- function (forward=TRUE,base.xvec,work.xvec,beta.vector,
ncats,nvars.total)
{
n.ones <- sum(base.xvec==1);  # indicates unique parameters
n.mults <- sum(base.xvec>1);  # indicates parameters constrained equal
if (forward) {
p <- 0;
if (n.ones>0) {
for (j in 1:ncats) {
for (i in 1:nvars.total) {
if (base.xvec[i,j] == 1) {
p <- p + 1;
beta.vector[p] <- work.xvec[i,j];
} #end of if
} #end of i
} #end of j
}
if (n.mults>0) {
nidxvals <- length(idxvals <- sort(unique(base.xvec[base.xvec>1])));
for (k in 1:nidxvals) {
for (j in 1:ncats) {
if (any(jktest <- base.xvec[,j]==idxvals[k])) {
beta.vector[p+k] <- work.xvec[jktest,j][1];
break;
}
}
}
}
return(beta.vector);
} else {
p <- 0;
if (n.ones>0) {
for (j in 1:ncats) {
for (i in 1:nvars.total) {
if (base.xvec[i,j] == 1) {
p <- p + 1;
work.xvec[i,j] <- beta.vector[p];
} #end of if
} #end of i
} #end of j
}
if (n.mults>0) {
nidxvals <- length(idxvals <- sort(unique(base.xvec[base.xvec>1])));
for (k in 1:nidxvals) {
for (j in 1:ncats) {
if (any(jktest <- base.xvec[,j]==idxvals[k])) {
work.xvec[jktest,j] <- beta.vector[p+k];
}
}
}
}
return(work.xvec);
} # end of else
} #end of mnl.xvec.mapping

############################################################################
## Create the jacstack (from tanh)                                     #
############################################################################

# jacstack.function:  arrange regressors for computing the Jacobian matrix
##   produces jacstack:  array of regressors,
##      dim(jacstack) = c(n observations, n UNIQUE parameters, n categories)
jacstack.function <- function(X,tvars.unique,xvec) {
xdim  <- dim(X)
obs   <- xdim[1]
#  nvars <- xdim[2]
ncats <- xdim[3]

jacstack <- array(0,dim=c(obs,tvars.unique,ncats));
n.ones <- sum(xvec==1);  # indicates unique parameters
n.mults <- sum(xvec>1);  # indicates parameters constrained equal
itu <- 0;
if (n.ones>0) {
for (j in 1:ncats) {
nxj <- sum(xvec[,j]==1);
if (nxj>0) {
jacstack[1:obs, itu + 1:nxj, j] <- X[rep(TRUE,obs),xvec[,j]==1,j];
}
itu <- itu + nxj;
}
}
if (n.mults>0) {
nxvecrows <- dim(xvec)[1];
nidxvals <- length(idxvals <- sort(unique(xvec[xvec>1])));
for (k in 1:nidxvals) {
for (j in 1:ncats) {
if (any(jktest <- xvec[,j]==idxvals[k])) {
kidx <- (1:nxvecrows)[jktest];  # xvec[,j] rows matching constraint k
nxj <- sum(jktest);
for (kk in 1:nxj) {
jacstack[1:obs, itu + k, j] <-
jacstack[1:obs, itu + k, j] + X[rep(TRUE,obs),kidx[kk],j];
}
}
}
}
}
return(jacstack)
} #jacstack.function

# jacstack.singles:  check jacstack array for regressors with a distinct value
#                    at only one observation
jacstack.singles <- function(jacstack) {
nunique <- dim(jacstack)[2];

jsingle <- rep(FALSE, nunique);
for (i in 1:nunique) {
if (length(unique(c(jacstack[,i,]))) == 2) {
jsingle[i] <- ifelse(any(table(c(jacstack[,i,])) == 1), TRUE, FALSE)
}
}
return(jsingle)
} #jacstack.singles.function

## 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.