[R] logistic regression and 3PL model

[John Fox]

Dear Mike,

You appear to have misinterpreted Brian Ripley's advice. The logitreg()
function that you've included in your message hasn't been modified to
specify a lower bound (higher than 0) to the probability of success. The
"lower" argument, which gets passed through to nlminb(), specifies the lower
bounds for the parameters, not for the fitted probabilities.

You also apparently have sent the S-PLUS version of the function. To get
this to run in R, you should replace the call to nlminb() with one to
optim(), and you'll have to modify the local function gmin() slightly. 

I hope this helps,
 John

--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
-------------------------------- 

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Michael Y. Lau
> Sent: Saturday, February 19, 2005 9:46 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] logistic regression and 3PL model
> 
> Hello colleagues,
> 
>  
> 
> This is a follow up to a question I posed in November 
> regarding an analysis I was working on.  Thank you to Dr. 
> Brian Ripley and Dr. John Fox for helping me out during that time.
> 
>  
> 
> I am conducting logistic regression on data set on psi (ESP) 
> ganzfeld trials.  The response variable is binary 
> (correct/incorrect), with a 25% guessing base rate.  Dr. 
> Ripley suggested that I modify a maximum likelihood fitting 
> of a binomial logistic regression presented in MASS$ (p. 445):
> 
>  
> 
> logitreg <- function(x, y, wt=rep(1, length(y)), intercept = 
> T, start=rep(0,p), ...) 
> 
> {
> 
>    fmin <- function (beta, X, y, w){
> 
>        p <- plogis(X %*% beta)
> 
>        -sum(2 * w * ifelse(y, log(p), log(1-p)))
> 
>        
> 
>    }
> 
>    gmin <- function (beta, X, y, w) {
> 
>        eta <- X %*% beta; p <- plogis(eta)
> 
>        -2 * (w *dlogis(eta) * ifelse(y, 1/p, -1/(1-p))) %*% X
> 
>    }
> 
>    if(is.null(dim(x))) dim(x) <- c(length(x), 1)
> 
>        dn <- dimnames(x)[[2]]
> 
>        if(!length(dn)) dn <- paste("Var", 1:ncol(x), sep="")
> 
>           p <- ncol(x) + intercept
> 
>           if(intercept) {x <- cbind(1,x); dn <- c("(Intercept)", dn)}
> 
>               if(is.factor(y)) y <- (unclass(y) !=1)
> 
>                  fit <- nlminb(start, fmin, gmin, X=x, y=y, 
> w=wt, method = "BFGS", ...)
> 
>                  names(fit$par) <- dn
> 
>                  cat("\nCoefficients:\n"); print(fit$par)
> 
>                  # R: use fit$value and fit$convergence
> 
>                  cat("\nResidual Deviance:", 
> format(fit$objective), "\n")
> 
>                  cat("\nConvergence message:", fit$message, "\n")
> 
>                  invisible(fit)
> 
> }
> 
>  
> 
> I have been using "lower=.25" in the call for logitreg:
> 
>  
> 
> options(contrasts = c("contr.treatment", "contr.poly"))
> 
> X <- model.matrix(longhit ~ age*gender, data=data)[,-1]
> 
> logitreg(X, data$longhit, lower =.25)
> 
> However, this is giving me .25 as the value for intercepts 
> and all the regression weights.  Am I correcting for the 
> guessing base rate correctly?  
> 
> Any suggestions or pointers would be greatly appreciated.  
> Thank you in advance.
> 
> Mike Lau
> 
>  
> 
> __________________________________
> Michael Y. Lau, M.A. 
> 118 Haggar Hall
> Department of Psychology
> University of Notre Dame
> Notre Dame, IN 46556
> (574) 631-1904
> mlau at nd.edu 
>   
> 
>  
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html