[R] Fitting a modified logistic with glm?

[Rubén Roa-Ureta]

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;                                # estimated 
proportion 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