Solved – Problems plotting GLM data of binomial proportional data

Question

This might be a question for the programmers, but I thought I would ask here first. Im comparing browsing pressure on plants between sites. Ive a value to indicate Browsing Pressure and count data of trees that have been damaged between locations. Ive been using Crawleys R example (page 574) regarding sex ratios and my issue is with the lines command to plot the information. I can't get a fitted regression line to the points to show up for my data. At least not using the commands Crawley suggested. Ive included Crawley's script for quick reference:

density = c(1,4,10,22,55,121,210,444)
females=c(1,3,7,18,22,41,52,79)
males=c(0,1,3,4,33,80,158,365)
numbers=as.data.frame(cbind(density,females,males))
attach(numbers)
par(mfrow=c(1,2))
p<-males/(males+females)
plot(log(density),p,ylab="Proportion male")
y<-cbind(males,females)
model<-glm(y~log(density),binomial)
xv<-seq(0,6,0.1)
plot(log(density),p,ylab="Proportion male")
lines(xv,predict(model,list(density=exp(xv)),type="response"))

My Script:

data<-####see the dput() data below to insert here
attach(data)
p1<-MS1/(MS1+M1)
y<-cbind(MS1,M1)
GM1<-glm(y~BPT,family=binomial (logit),data=data)
summary(GM1)
plot(BPT,p1,col="black",pch=1,main="Relationship a",xlab="Browsing pressure",
ylab="Moose Damage Survey")

This is where I am having difficulties interpreting Crawley's lines command. I can create a sequence for my data like his xv values – mine go from 0 to 0.7.
xv<-seq(0,0.7,0.01) But as Im not sure what is happening in the command line, blindly inserting my data results in no line showing on my graph. I have some success with regLine from the car library, but its not a line fit to my data. Any help explaining why this works for Crawleys but not my data would be appreciated

My DATA:

structure(list(S = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 
30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 
46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60), 
    BPT = c(0.195778643105884, 0.0651216427326909, 0.0199432997190648, 
    0.255717971810655, 0.21031730503132, 0.0380060076389563, 
    0.0592940237882193, 0.00425756881206537, 0.107351336677244, 
    0.353626952501855, 0.134358354834345, 0.0774003112699534, 
    0.0571275803730281, 0.102480049705623, 0.506048974065245, 
    0.637795327349053, 0, 0.00825369912299262, 0.123063816485164, 
    0.244214078077056, 0, 0.231982659809539, 0.0265824936284919, 
    0, 0, 0, 0.197470272555566, 0, 0.0481572367104686, 0.339072491473947, 
    0.0640004726042278, 0.099465609734403, 0, 0.153951519788959, 
    0.103791449528225, 0.22066571012681, 0.101286659375451, 0.0210055739071239, 
    0, 0, 0.00811999120846242, 0.0576334384710714, 0.0514027186783361, 
    0.00907843357263548, 0.0785740169740806, 0.0801612249152947, 
    0.0439884643248279, 0.189608233745472, 0.434349352638054, 
    0.188190501082836, 0.016381369054353, 0.215189432724754, 
    0.0948563822493123, 0, 0.0111745795351048, 0.0289751701436048, 
    0.0962515835377456, 0.199970627051051, 0.0306277369169669, 
    0), M1 = c(51, 123, 137, 23, 69, 98, 80, 59, 84, 87, 63, 
    88, 64, 75, 30, 19, 86, 152, 55, 22, 115, 66, 75, 124, 87, 
    179, 39, 66, 117, 37, 59, 103, 57, 65, 83, 68, 87, 104, 134, 
    89, 52, 61, 65, 75, 71, 97, 45, 37, 86, 82, 36, 139, 40, 
    56, 63, 64, 37, 79, 70, 121), M2 = c(108, 195, 255, 31, 104, 
    157, 135, 88, 117, 152, 96, 124, 102, 118, 37, 25, 135, 277, 
    97, 47, 174, 121, 134, 196, 163, 270, 75, 123, 214, 46, 109, 
    156, 92, 118, 127, 112, 134, 152, 194, 132, 107, 101, 105, 
    140, 107, 163, 87, 60, 122, 122, 66, 237, 66, 95, 108, 95, 
    50, 128, 128, 189)), .Names = c("S", "BPT", "M1", "M2"), row.names = c(NA, 
-60L), class = "data.frame")

Answer 1

Here is a solution.

First, I replaced MS1 with M2 since there is no varibale with the name MS1 in your example data.

attach(data)
p1<-M2/(M2+M1)
y<-cbind(M2,M1)

The model:

GM1<-glm(y~BPT,family=binomial (logit),data=data)

Note that there is one important difference between your model and Crawley's. His predictor varibale was log-transformed, your predictor variable BPT is used untransformed. By the way: Since BPT does include zeros, log-transformation would result in -Inf values.

Now the main plot:

plot(BPT,p1,col="black",pch=1,main="Relationship a",xlab="Browsing pressure",
     ylab="Moose Damage Survey")

Now the important differences. Crawley created a sequence xv containg values corresponding to his log-transformed predictor variable density. Note that the lengths of density and xv must be identical. Compare:

range(log(density))
#[1] 0.000000 6.095825
xv<-seq(0,6,0.1)
range(xv)
#[1] 0 6

Since your predictor variable is untransformed, you don't need this kind of sequence but can use the original variable BPT.

lines(BPT,predict(GM1,type="response"))

This is the plot:

If you would like to use your own sequence of values for the prediction, use this:

xv<-seq(0,0.7,length.out = length(BPT))
lines(xv,predict(GM1, list(BPT = xv), type="response"))