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@(protected)
[[alternative HTML version deleted]]
______________________________________________
R-help@(protected)
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
