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
[EMAIL PROTECTED]
http://haddock.bordeaux.cemagref.fr:8080/ocms1/opencms/estuaires/Poissons_migrateurs_et_changement_global.html



        [[alternative HTML version deleted]]

______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to