Dear all,
To modelize the abundance of fish (4 classes) with a set of environmental
variables, I used the polr and predict.polr functions. I would like to know how
to bring the cumulated probabilities back to a discrete ordinal scale.
For the moment I used the predict.polr function with the argument
"class". Is there an other way?
polrf <- polrf <- polr_mod(formula = acipenser_gueldenstaedtii ~ Long +
poly(Surf, 2, raw = TRUE) + poly(TempSum, 2, raw = TRUE) , data = mydata, method
= "logistic", Hess = TRUE, na.action = na.omit)
pred1 <- predict.polr(polrf, newdata = mydata2, type = "class")
predict.polr <- function(object, newdata,
type=c("class","probs"), ...)
{
if(!inherits(object, "polr")) stop("Not a polr fit")
type <- match.arg(type)
if(missing(newdata)) Y <- object$fitted
else {
newdata <- as.data.frame(newdata)
Terms <- delete.response(object$terms)
m <- model.frame(Terms, newdata, na.action = function(x) x)
X <- model.matrix(Terms, m, contrasts = object$contrasts)
xint <- match("(Intercept)", dimnames(X)[[2]], nomatch=0)
if(xint > 0) X <- X[, -xint, drop=F]
n <- nrow(X)
q <- length(object$zeta)
eta <- drop(X %*% object$coef)
cumpr <- matrix(plogis(matrix(object$zeta, n, q, byrow=T) - eta), ,
q)
Y <- t(apply(cumpr, 1, function(x) diff(c(0, x, 1))))
dimnames(Y) <- list(dimnames(X)[[1]], object$lev)
}
if(!is.null(object$na.action)) Y <- napredict(object$na.action, Y)
switch(type, class={
Y <- factor(max.col(Y), levels=seq(along=object$lev),
labels=object$lev)
}, probs={})
drop(Y)
}
Thank you.
Géraldine LASSALLE
Cemagref
Unité Ecosystèmes estuariens et
Poissons migrateurs amphihalins
50 avenue de Verdun
33612 Gazinet-Cestas
Tel: 05.57.89.09.98
Fax: 05.57.89.08.01
geraldine.lassalle@bordeaux.cemagref.fr
http://haddock.bordeaux.cemagref.fr:8080/ocms1/opencms/estuaires/Poissons_migrateurs_et_changement_global.html
[[alternative HTML version deleted]]