[19155] in s-news-athena
Re: [S] Conditional Multinomial Logit
daemon@ATHENA.MIT.EDU (blparis@kodak.com)
Wed Sep 8 17:44:44 1999
From: blparis@kodak.com
To: s-news@wubios.wustl.edu
Message-Id: <852567E6.00761DDE.00@knotes.kodak.com>
Date: Wed, 8 Sep 1999 17:40:02 -0400
Mime-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline
From: Brandon L. Paris
Thanks to all who replied to my original post about conditional logit model
in S-plus. What follows is a summary of the two replies. It turns out
that both require the same data manipulation I am presently doing. Michael
Conklin provides code that requires that the data already be structured as
1 row per product (per Choice set, I presume). Brad Biggerstaff provides
wrapper code to appropriately format the data before calling coxph(). I
haven't had the chance to try them, but will attempt to do so in the next
few days.
I have just one comment/question about the submissions. I have read (in
several sources) that the only appropriate method for fitting this type of
model is to use the "Breslow" method. Brad's code uses the "Exact" method.
Does anyone have any info about which is more appropriate for Discrete
Choice models? Thanks.
Michael Conklin wrote:
The data set is set up with 1 row per product. The Choice variable is
coded 1 if the product was chosen and 0 if it was not. The Choice set
number is also in the data and the appropriate product descriptor
variables. This leads to a data set with (Num Choice Sets per
Respondent)* (Num Products per choice set) * Num Respondents rows.
Then I use the following code.
survchoice<-Surv(2-choicedataset$Choice,is.element(choicedataset$Choice,1))
options(contrasts=c("contr.treatment","contr.treatment"))
totalmodel<-coxph(survchoice~V1*V2+V3+V4+V5+V6+strata(ChoiceSet),data=choic
edataset,method="breslow")
summary(totalmodel)
***************************************************************************
***************************
Brad Biggerstaff wrote:
I don't know about your specific application, but with regards to the
manipulations required to use coxph, I wrote a "wrapper function" to do
matched logistic regression (conditional logistic) which takes as its input
the standard glm-type arguments plus a parameter for the matching variable.
The returned object is of class "coxph" so the methods for it work for
this...e.g., print. Here it is, and you can scan the code to find the
manipulations you're after; they're, in particular, at the lines where
PseudoTime and PseudoSurv are set and that the call to coxph uses
method="exact".
Cheers,
Brad
Brad Biggerstaff, Ph.D.
Centers for Disease Control and Prevention
National Center for Infectious Diseases
Division of Vector-Borne Infectious Diseases
P.O. Box 2087
Fort Collins, CO 80522-2087 U.S.A.
(970) 221-6473 bkb5@cdc.gov
"mlogitreg"<-
function(formula = formula(data), data = sys.parent(), weights, subset = T,
na.action, matching.variable = NULL, eps = 0.0001, iter.max = 10, x
= F,
y = T)
{
if(is.null(matching.variable))
stop("matching.variable must be specified")
call <- match.call()
m <- match.call()
m$matching.variable <- m$eps <- m$iter.max <- m$x <- m$y <- NULL
# m$dataset <- as.character(unlist(m)["data"])
m[[1]] <- as.name("model.frame")
m <- eval(m, sys.parent())
Terms <- attr(m, "terms")
CaseStatus <- model.extract(m, "response")
weights <- model.extract(m, "weights")
offset <- attr(Terms, "offset")
if(length(unique(CaseStatus)) != 2)
stop(" Response must have only two levels")
if(is.factor(CaseStatus)) {
CaseStatus <- codes(CaseStatus)
CaseStatus[CaseStatus == 2] <- 0
}
if(options()$contrasts[1] != "contr.treatment") {
options(contrasts = c("contr.treatment", "contr.poly"))
cat("\n ...factor contrasts reset to contr.treatment\n\n")
}
X <- model.matrix(Terms, m)
PseudoTime <- rep(1, length(CaseStatus))
PseudoSurv <- Surv(PseudoTime, CaseStatus)
if(is.data.frame(data)) {
if(!is.character(matching.variable[[1]]))
MatchingVar <- data[subset, as.character(
matching.variable[[1]])]
else MatchingVar <- data[subset, matching.variable[[1]]]
}
else MatchingVar <- eval(as.name(matching.variable[[1]]))[subset]
fit <- coxph.fit(X, PseudoSurv, strata = strata(MatchingVar), offset
=
offset, weights = weights, method = "exact", eps = eps,
iter.max = iter.max, rownames = row.names(m))
fit$coefficients <- fit$coefficients[-1] # get rid of
Intercept
fit$var <- fit$var[-1, -1] #get rid of Intercept row and column
attr(fit, "class") <- "coxph"
fit$n <- length(CaseStatus)
fit$call <- call
if(x)
fit$x <- X
if(y)
fit$y <- CaseStatus
fit
}
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
