Hey all.
Question category: R and regression; Posted in R forum since R code.
Question: How best to quantify variation in an inflection point for a logistic curve?
Background of problem: Birdies begin their birdie lives contained in eggs. Then, they hatch and continue their lives. Bothersome scientists are interested in determining how many days the birdies spend as eggs; that is, how many days until they hatch. We check nests each 2-4 days and record 1=eggs present, 2=hatching. Thus, because we often don't catch the birdies as they're emerging from their shells, there is some uncertainty in observed hatch date.
Current modeling approach:(in line with the way I already have my data for other analyses) Binomial GLM: x=nest age (in days), response: 1=eggs, 0=nestlings. Model the heck out of it.
Generate example and plot:
# Example plot
Estimation: From model equation, use beta and SE.beta to back-calculate x at y=0.5 (point of inflection): logit(0.5)=intercept+beta*x. Draw n=10000 from normal distribution on link scale.
I've tried other approaches, but I think this may do a good job of both giving me a good estimate of hatch day while also allowing for observer uncertainty.
Thoughts? Anything I'm not considering or approached with less-then-impressive rigor? (I always worry about drawing from the model output distributions, but it always seem to come out close in tests.)
Question category: R and regression; Posted in R forum since R code.
Question: How best to quantify variation in an inflection point for a logistic curve?
Background of problem: Birdies begin their birdie lives contained in eggs. Then, they hatch and continue their lives. Bothersome scientists are interested in determining how many days the birdies spend as eggs; that is, how many days until they hatch. We check nests each 2-4 days and record 1=eggs present, 2=hatching. Thus, because we often don't catch the birdies as they're emerging from their shells, there is some uncertainty in observed hatch date.
Current modeling approach:(in line with the way I already have my data for other analyses) Binomial GLM: x=nest age (in days), response: 1=eggs, 0=nestlings. Model the heck out of it.
Generate example and plot:
Code:
# generate data
require(boot)
x<-1:100
set.seed=(1001)
mid.vals<-rbinom(20,size=1,prob=0.5)
y<-c(rep(1,40), mid.vals, rep(0,40))
mod<-glm(y~x, family=binomial)
pred.link<-predict(mod,se.fit=T)$fit
pred.se<-predict(mod,se.fit=T)$se.fit
pred.val<-inv.logit(pred.link)
pred.upper<-inv.logit(pred.link+1.96*pred.se)
pred.lower<-inv.logit(pred.link-1.96*pred.se)
Code:
plot(x,pred.val, type="l",lwd=2, col="white",xlab="",ylab="",ylim=c(0,1))
polygon(c(x,rev(x)),c(pred.lower,rev(pred.upper)),col="grey80",border=NA)
par(new=T)
plot(x,pred.val, type="l",lwd=2,xlab="",ylab="",ylim=c(0,1))
par(new=T)
plot(x,pred.upper, type="l",lwd=2,lty=2,xlab="",ylab="",ylim=c(0,1))
par(new=T)
plot(x,pred.lower, type="l",lwd=2,lty=2,xlab="nest age (days)",ylab="Pr(unhatched egg)",ylim=c(0,1))
abline(h=0.5, col="red",lwd=1) # horizontal line at y=0.5
abline(v=49.5, col="red",lwd=1) # general line to see x-point of inflection
Code:
inter.draws<-rnorm(10000, coefficients(mod)[1],sqrt(diag(vcov(mod)))[1])
slope.draws<-rnorm(10000, coefficients(mod)[2],sqrt(diag(vcov(mod)))[2])
mean.hatch.day<-inter.draws/(-1*slope.draws) #to estimate mean x from y=0.5, or logit(0.5)
hist(mean.hatch.day)
quantile(mean.hatch.day)
Thoughts? Anything I'm not considering or approached with less-then-impressive rigor? (I always worry about drawing from the model output distributions, but it always seem to come out close in tests.)