Solved – Problem plotting GLM data of binomial proportional data 2

binomial distributiondata visualizationgeneralized linear modelproportion;r

I am trying to plot the predicting line for a GLM with proportional data. After having troubles with the example from "Crawley" I found this link:
Problems plotting GLM data of binomial proportional data
My data consists out of 3 columns: Altitude, number of beetles and number of other insects. I want to investigate the change of beetle proportion throughout the different altitudes.
After importing the table I calculated the proportion as shown in the topic linked above.
Then

y=cbind(mean$beetles,mean$rest)   
model=glm(y~altitude,family=binomial(logit),data=data)   
plot(p~data$altitude, ylab="Beetles (proportional)", xlab="Altitude")
p is the propotion i calculated before
lines(mean$altitude,predict(model,type="response"))

I then get this plot:

plot

this is obviously not correct but i cannot figure out what is wrong here.
I would be very pleased if you helped me.

this is the little data set:

 altitude,beetles,rest
 304,88,35
 842,7,18
 877,21,27
 1505,2,25
 1118,11,19
 702,25,24
 806,16,17
 960,16,32
 975,34,61
 1112,11,17
 1348,11,32
 502,58,24
 1428,2,13
 1521,0,35
 1040,9,6
 1412,8,19
 1231,5,19
 980,14,16
 1053,2,10

Best Answer

Your question does not specify how the p was calculated, but I assume it is the percentage of beetles on every plot you have sampled. I can't see the image you've attached either (I'm probably out of luck for today), but judging from the code you've posted, I'd guess you need to sort the predicted values before plotting them with lines(). Something like this should work:

# This is the data, saved in object d
d<-structure(list(altitude = c(304L, 842L, 877L, 1505L, 1118L, 702L, 
806L, 960L, 975L, 1112L, 1348L, 502L, 1428L, 1521L, 1040L, 1412L, 
1231L, 980L, 1053L), beetles = c(88L, 7L, 21L, 2L, 11L, 25L, 
16L, 16L, 34L, 11L, 11L, 58L, 2L, 0L, 9L, 8L, 5L, 14L, 2L), rest = c(35L, 
18L, 27L, 25L, 19L, 24L, 17L, 32L, 61L, 17L, 32L, 24L, 13L, 35L, 
6L, 19L, 19L, 16L, 10L)), .Names = c("altitude", "beetles", "rest"
), class = "data.frame", row.names = c(NA, -19L))

# This fits the model
model=glm(cbind(d$beetles, d$rest)~altitude,family=binomial(logit),data=d)   

# The plot of proportion of beetles against the altitude
plot(prop.table(cbind(d$beetles, d$rest), 1)[,1]~d$altitude, ylab="Beetles (proportional)", xlab="Altitude")

# The prediction of the proportion of the beetles from the fit
pred<-predict(model,type="response")

# Ordering the predicted values (and altitude) according to the altitude
pred2<-pred[order(d$altitude)]
    alt<-d$altitude[order(d$altitude)]

# Drawing the line
lines(alt,pred2)

Or, a simpler method, sort the data frame first:

d<-d[order(d$altitude),]
model=glm(cbind(d$beetles, d$rest)~altitude,family=binomial(logit),data=d)   
plot(prop.table(cbind(d$beetles, d$rest), 1)[,1]~d$altitude, ylab="Beetles (proportional)", xlab="Altitude")
lines(d$altitude,predict(model,type="response"))
Related Question