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") Best Answer
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"))