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

Defines functions mnl.probfuncres.stdkth.smallest

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

# probfunc: matrix of estimated probabilities
mnl.probfunc <-  function(Y, Ypos, Xarray, tvec) {
nobs <- dim(Y)[1]
ncats <- dim(Y)[2]
eta <- matrix(0,nobs,ncats)
for (j in 1:ncats) {
useobs <- Ypos[,j];
if (dim(tvec)[1] == 1) {
eta[useobs,j] <- exp(Xarray[useobs,,j] * tvec[,j]);
}
else {
eta[useobs,j] <- exp(Xarray[useobs,,j] %*% tvec[,j]);
}
}
return( c(1/(eta %*% rep(1,ncats))) * eta )
}#end of mnl.probfunc

## residual.generator:  raw and median-centered ortho-standardized residuals
residual.generator <- function (tvec, Y, Ypos, X, m)
{

y.prob <- mnl.probfunc(Y, Ypos, X, tvec);

if (all(Ypos)) {
Sres.raw <- res.std(Y, m, y.prob);
Sres <- Sres.raw - median(Sres.raw);
}
else {
nobs <- dim(Y)[1];
ncats <- dim(Y)[2];
Sres.raw <- matrix(0, nobs, ncats-1);
hasall <- apply(Ypos, 1, sum) == ncats;
nobsall <- sum(hasall);
if (nobsall > 0) {
Yuse <- matrix(Y[hasall,], nobsall, ncats);  # in case nobsall == 1
puse <- matrix(y.prob[hasall,], nobsall, ncats);
Sres.raw[hasall,] <- res.std(Yuse, c(Yuse %*% rep(1,ncats)), puse);
}
hasless <- (1:nobs)[!hasall];
Sres.use <- matrix(TRUE, nobs, ncats-1);
for (i in hasless) {
usecats <- Ypos[i,];
nlesscats <- sum(usecats);
ocats <- 1:(nlesscats-1);
Yuse <- matrix(Y[i,usecats], 1, nlesscats);
puse <- matrix(y.prob[i,usecats], 1, nlesscats);
Sres.raw[i,ocats] <- res.std(Yuse, c(Yuse %*% rep(1,nlesscats)), puse);
Sres.use[i,nlesscats:(ncats-1)] <- FALSE;
}
Sres <- Sres.raw - median(Sres.raw[Sres.use]);
Sres[!Sres.use] <- 0;
}

return(list(Sres.raw=Sres.raw,Sres=Sres));
}

## res.std: compute orthogonalized and standardized residuals
##  based on Tanabe and Sagae (1992 JRSSB)
res.std <- function(y,m,p, print.level=0)
{
obs    <- dim(y)[1]
ncats  <- dim(y)[2]
summat <- matrix(0,ncats,ncats)
for (i in 1:(ncats-1)) summat[(i+1):ncats,i] <- 1;
q <- p %*% summat
## d:  nonzero values in the diagonal matrix of the Cholesky decomposition
d <- matrix(0,obs,ncats-1)
for (i in 1:(ncats-1)) {
if (i==1) d[,i] <- p[,i]*q[,i];
if (i>1)  d[,i] <- p[,i]*q[,i]/q[,i-1];
}
## r: raw residuals
r <- y-m*p;
summat <- matrix(0,ncats,ncats)
for (i in 1:ncats) summat[i,i:ncats] <- 1;
rsum <- r %*% summat
## Or:  orthogonalized residuals
Or <- matrix(0,obs,ncats-1)
for (i in 1:(ncats-1)) {
if (i==1) Or[,i] <- r[,i];
if (i>1)  Or[,i] <- r[,i] + rsum[,i-1]*p[,i]/q[,i-1];
}
## Sr: standardized residuals
Sr <- Or/sqrt(m*d);

if (print.level > 0) {
cat("res.std: Or\n");
print(Or);
cat("res.std: m\n");
print(m);
cat("res.std: d\n");
print(d);
}
return(Sr);
}

kth.smallest  <- function(SortVector, obs, k)
{
if (obs < k) {
print("ERROR! obs < k in kth.smallest\n");
return(-1);
}
return(.Call("kthSmallest",
as.double(SortVector), as.integer(obs), as.integer(k),
PACKAGE="multinomRob"));
} #end of kth.smallest
```

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.