############################################
#### Segmented Regression Demo ####
############################################

#install.packages("MASS")
#install.packages("segmented")
library(MASS) #for dose.p function
library(segmented)

#Predictor variable 'V'
V = c(0,25,45,20,37,30,10,35,10,60,55,90,80,45,85,63,60,75,53,100)
#Present / Absent response variable 'P_A'
P_A = c(0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1)
#plot binomial data
plot(P_A~V)

#logistic regression analysis of raw data
P_A.V.glm <- glm(P_A~V,  family=binomial)
summary(P_A.V.glm)#significant difference between P and A
dose.p(P_A.V.glm, p=0.5) #quantifies 'V' the 'midpoint'

#predicted data for the logistic curve
range(V)
V.pred <- seq(0,100,1)
P_A.pred <- predict(P_A.V.glm, list(V =V.pred), type="response")
#Plot logistic curve
plot(P_A.pred~V.pred) #nice lazy 'S' shaped curve

#Generate a model of the logistic curve for 'segmented' to work from
#(Here we are analyzing the shape of the curve, not the raw data)
P_A.pred.mod <- lm(P_A.pred~V.pred) 
#Run the segmented regression
P_A.pred.segmented <- segmented(P_A.pred.mod, seg.Z = ~V.pred, psi=c(10,60))  
#psi sets the start of the section(s) of the curve that you want to analyse. 
summary(P_A.pred.segmented) #identifies the inflection points  

#Breakpoint plot
plot(P_A.pred~V.pred,xlab="Predictor Variable 'V'", ylab="Response Variable 'P_A'", cex=0,
     main="Lower(Concave) and Upper(Convex) Inflection Points on the Logistic Curve")
axis(side=2, at=seq(0,1.0, by=0.1))
lines(V.pred, P_A.pred, col="grey", lwd=3)
plot(P_A.pred.segmented, add=T, col="blue", lwd=3)
#annotations and line segments
segments(x0=62.259, x1=62.259, y0=0, y1=0.9, lwd=1, lty=11)
segments(x0=0, x1=62.259, y0=0.9, y1=0.9, lwd=1, lty=11)
text(x=65.3,y=0.2, label="Breakpoint = 62.3", srt=-90)
segments(x0=48.836, x1=48.836, y0=0, y1=0.5, lwd=1, lty=11, col="grey")
segments(x0=0, x1=48.836, y0=0.5, y1=0.5, lwd=1, lty=11, col="grey")
text(x=52,y=0.2, label="Midpoint = 48.8", srt=-90)

########################################################################