From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>

Date: Sat 16 Jul 2005 - 01:58:51 EST

Date: Sat 16 Jul 2005 - 01:58:51 EST

Robin Hankin <r.hankin@noc.soton.ac.uk> writes:

*> Hi
**>
*

> I am trying to make glm() work to analyze a toy logit system.

*>
**> I have a dataframe with x and y independent variables. I have
**>
**> L=1+x-y (ie coefficients 1,1,-1)
**>
**> then if I have a logit relation with L=log(p/(1-p)),
**> p=1/(1+exp(L)).
**>
**> If I interpret "p" as the probability of success in a Bernouilli
**> trial, and I can observe the result (0 for "no", 1 for "yes")
**> how do I retrieve the coefficients c(1,1,-1)
**> from the data?
**>
**> n <- 300
**> des <- data.frame(x=(1:n)/n,y=sample(n)/n) # experimental design
**> des <- cbind(des,L=1+des$x-des$y) # L=1+x-y
**> des <- cbind(des,p=1/(1+exp(des$L))) # p=1/(1+e^L)
**> des <- cbind(des,obs=rbinom(n,1,des$p)) # observation: prob of
**> success = p.
**>
**>
**> My attempt is:
**>
**> glm(obs~x+y,data=des,family=binomial(link="logit"))
**>
**> But it does not retrieve the correct coefficients of c(1,1,-1) ;
**> I would expect a reasonably close answer with so much data.
**>
**> What is the correct glm() call to perform my logit analysis?
*

Apart from a sign error, the only fault seems to be that you think that 300 is a large number... Try upping n to 30000 and you'll see.

(The sign error is of course that log(p/(1-p)) = L implies that p = exp(L)/(1+exp(L)) = 1/(1+exp(-L))).

-- O__ ---- Peter Dalgaard ุster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907 ______________________________________________ R-help@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.htmlReceived on Sat Jul 16 02:02:20 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:33:45 EST
*