From: Rubén Roa-Ureta <rroa_at_udec.cl>

Date: Sat, 08 Nov 2008 16:59:45 -0300

negloglik <- -sum(count*(mat*log(prop.est)+imat*log(iprop.est))); #binomial likelihood, prob. model

}

prop.lik <- nlm(fn,p=c(25.8,33.9),hessian=TRUE)

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sat 08 Nov 2008 - 20:03:08 GMT

Date: Sat, 08 Nov 2008 16:59:45 -0300

Mike Lawrence wrote:

> Hi all,

*>
**> Where f(x) is a logistic function, I have data that follow:
**> g(x) = f(x)*.5 + .5
**>
**> How would you suggest I modify the standard glm(..., family='binomial')
**> function to fit this? Here's an example of a clearly ill-advised attempt to
**> simply use the standard glm(..., family='binomial') approach:
**>
**> ########
**> # First generate some data
**> ########
**> #define the scale and location of the modified logistic to be fit
**> location = 20
**> scale = 30
**>
**> # choose some x values
**> x = runif(200,-200,200)
**>
**> # generate some random noise to add to x in order to
**> # simulate real-word measurement and avoid perfect fits
**> x.noise = runif(length(x),-10,10)
**>
**> # define the probability of success for each x given the modified logistic
**> prob.success = plogis(x+x.noise,location,scale)*.5 + .5
**>
**> # obtain y, the observed success/failure at each x
**> y = rep(NA,length(x))
**> for(i in 1:length(x)){
**> y[i] = sample(
**> x = c(1,0)
**> , size = 1
**> , prob = c(prob.success[i], 1-prob.success[i])
**> )
**> }
**>
**> #show the data and the source modified logistic
**> plot(x,y)
**> curve(plogis(x,location,scale)*.5 + .5,add=T)
**>
**> ########
**> # Now try to fit the data
**> ########
**>
**> #fit using standard glm and plot the result
**> fit = glm(y~x,family='binomial')
**> curve(plogis(fit$coefficients[1]+x*fit$coefficients[2])*.5+.5,add=T,col='red',lty=2)
**>
**> # It's clear that it's inappropriate to use the standard
**> "glm(y~x,family='binomial')"
**> # method to fit the modified logistic data, but what is the alternative?
**>
*

A custom-made fit function using nlm?

Below, in the second line the logistic model has been re-parameterized
to include p[1]=level of the predictor which predicts a 50% probability
of success, and p[2]=level of the predictor which predicts a 95%
probability of success. You can change the model in this line to your
version. In the 4th line (negative log-likelihood) mat is a vector of 0s
and 1s representing success and imat is a vector of 1s and 0s
representing failure. The fit is for grouped data.
fn <- function(p){

prop.est <- 1/(1+exp(log(1/19)*(l-p[1])/(p[2]-p[1]))); # estimated
proportion of success

iprop.est <- 1-prop.est; # estimatedproportion of failure

negloglik <- -sum(count*(mat*log(prop.est)+imat*log(iprop.est))); #binomial likelihood, prob. model

}

prop.lik <- nlm(fn,p=c(25.8,33.9),hessian=TRUE)

**HTH
**

Rubén

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sat 08 Nov 2008 - 20:03:08 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Sat 08 Nov 2008 - 21:30:24 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*