On Fri, 15 Jul 2005, Robin Hankin wrote: > 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)).
Not quite, see below. > 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) A logit would be p = e^L/(1+e^L), so your signs for L are reversed. > 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. You actually have so little data. > What is the correct glm() call to perform my logit analysis? The call is correct, the expectation is not. A single bernoulli observation provides far less information than you seem to suppose. I got Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.4747 0.3670 -4.019 5.85e-05 *** x -0.5549 0.4672 -1.188 0.23494 y 1.2963 0.4731 2.740 0.00614 ** and note how large the standard errors are. With 10000 examples you will get closer. Having fixed your sign change, I got Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.98711 0.06024 16.39 <2e-16 *** x 1.00896 0.08052 12.53 <2e-16 *** y -0.87798 0.08031 -10.93 <2e-16 *** -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595 ______________________________________________ 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.html