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?
--
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University
www.thatmike.com
Looking to arrange a meeting? Check my public calendar and find a time that
is mutually free:
http://tinyurl.com/mikescalendar
(tiny url redirects to google calendar)
~ Certainty is folly... I think. ~
[[alternative HTML version deleted]]
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
Mike Lawrence <mike <at> thatmike.com> writes:> 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: > #define the scale and location of the modified logistic to be fit > location = 20 > scale = 30 > x = runif(200,-200,200) > x.noise = runif(length(x),-10,10) > prob.success = plogis(x+x.noise,location,scale)*.5 + .5 > 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]) > ) > } > plot(x,y) > curve(plogis(x,location,scale)*.5 + .5,add=T)Hi, You might try the link mafc.logit(m = 2) defined in the psyphy package. Continuing with your example, library(psyphy) fit <- glm(y ~ x, binomial(mafc.logit(2)), control = glm.control(maxit = 100)) # default didn't converge x.ord <- order(x) lines(x[x.ord], fitted(fit)[x.ord], col = "red", lwd = 3) HTH, Ken -- Ken Knoblauch Inserm U846 Institut Cellule Souche et Cerveau D?partement Neurosciences Int?gratives 18 avenue du Doyen L?pine 69500 Bron France tel: +33 (0)4 72 91 34 77 fax: +33 (0)4 72 91 34 61 portable: +33 (0)6 84 10 64 10 http://www.sbri.fr
Hi Mike, the model you consider is a special case of the four-parameter logistic model where the lower and upper asymptotes are fixed at 0.5 and 1, respectively. Therefore, this (dose-response) model can fitted using the R package 'drc': library(drc) xy.m <- drm(y~x, fct = L.4(fixed=c(NA,0.5, 1, NA)), type = "binomial") plot(x,y) lines(x[x.ord], fitted(fit)[x.ord], col = "red", lwd = 3) # psyphy fit plot(xy.m, add = TRUE, log = "", col = 3, type = "fit", lwd = 3) Christian