Chad-
The plot you hope to create is meaningful when the variable temperature is
treated as a continuous variable. Below is some code that treats
temperature as a continuous variable and creates the plot. Note that
temperature enters the model "linearly", and you may want to explore other
possibilities for temperature beyond this one.
-tgs
Z <- textConnection("
Temp Sal Go Total
5 34 1 1
5 34 1 1
5 34 1 1
5 21 1 1
5 21 1 1
5 21 0 1
10 34 1 1
10 34 0 1
10 34 1 1
10 21 1 1
10 21 0 1
10 21 0 1
15 34 0 1
15 34 0 1
15 34 0 1
15 21 0 1
15 21 0 1
15 21 0 1")
ddd <- read.table(Z,header=TRUE)
close(Z)
#Fit model
glm1 <- glm(cbind(Go,Total-Go)~Temp*factor(Sal),data=ddd,family="binomial")
#Calculate predicted values
newdata1 <- data.frame(Temp = seq(5,15,length=50), Sal = rep(34,50))
pred34 <- predict(glm1,newdata=newdata1,type="response")
newdata2 <- data.frame(Temp = seq(5,15,length=50), Sal = rep(21,50))
pred21 <- predict(glm1,newdata=newdata2, type="response")
#Plot model predicted curves
plot(c(5,15),c(0,1),type="n",xlab="Temperature",ylab="Probability of going")
lines(newdata1$Temp,pred34,lwd=3,col="red")
lines(newdata2$Temp,pred21,lwd=3,col="blue")
#This next bit of code calculates the "observed" probabilities. I admit,
it isn't the clearest code.
ttt <- as.data.frame(xtabs(~Temp+Sal+Go,data=ddd)[,,2] /
xtabs(~Temp+Sal+Total,data=ddd)[,,1])
ttt$X <- as.numeric(rownames(ttt))
#Plot observed probabilities
points(ttt$X,ttt$"21",cex=4,lwd=2,pch=16,col="white")
points(ttt$X,ttt$"34",cex=4,lwd=2,pch=16,col="white")
points(ttt$X,ttt$"21",cex=4,lwd=2,col="blue")
points(ttt$X,ttt$"34",cex=4,lwd=2,col="red")
text(ttt$X,ttt$"21","21",col="blue")
text(ttt$X,ttt$"34","34",col="red")
On Wed, Dec 12, 2012 at 7:50 AM, Chad Widmer <cw...@st-andrews.ac.uk> wrote:
> Dear R Wizards,
>
> After much frustration and days of confusion I have finally broken down and
> am asking for help, which I dont like doing, but I just cant figure this
> one out on my own. Ive conducted a laboratory experiment testing the
> effects of temperature and salinity on whether or not a biological event
> will occur (Go or NoGo). Ive coded the factors temperature and salinity
> as factors for the binomial glm, and I havent had any trouble fitting the
> model and checking assumptions.
>
> I am however having trouble with the predict.glm function. I want to
> create a graph using my data that is similar to the one produced by the
> budworm example at the bottom of the predict.glm R documentation. In my
> case I want temperature on the x axis, probability on the y axis, and the
> lines on the graph to represent the probability of the event occurring at
> the different salinities tested at the different temperatures.
>
> I created a smaller version of my data and have included it and the R code
> I used below. I get two main problems that I hope youre willing to help
> with.
>
> 1. When I input the text argument the output in the graph gives salinity
> values that are superimposed on one another. And, the values dont seem to
> make sense they are returned as probabilities of either 1.0 or 0.
>
> 2. When I input the lines argument I get the following error messages:
>
> Error: variable 'fTemp' was fitted with type "factor" but type "numeric"
> was supplied
>
> In addition: Warning message:
>
> In model.frame.default(Terms, newdata, na.action = na.action, xlev =
> object$xlevels) : variable 'fTemp' is not a factor
>
> grrrrrrrrrr
>
> Pleasehelp<-read.table("Rhelp.txt",h=T)
>
> attach(Pleasehelp)
>
> fix(Pleasehelp)
>
> Temp Sal Go Total
>
> 5 34 1 1
>
> 5 34 1 1
>
> 5 34 1 1
>
> 5 21 1 1
>
> 5 21 1 1
>
> 5 21 0 1
>
> 10 34 1 1
>
> 10 34 0 1
>
> 10 34 1 1
>
> 10 21 1 1
>
> 10 21 0 1
>
> 10 21 0 1
>
> 15 34 0 1
>
> 15 34 0 1
>
> 15 34 0 1
>
> 15 21 0 1
>
> 15 21 0 1
>
> 15 21 0 1
>
> fTemp<-factor(Temp)
>
> fSal<-factor(Sal)
>
> Go<-Go
>
> NoGo<-Total-Go
>
> Went<-cbind(Go,NoGo)
>
> DF<-data.frame(fTemp,fSal,Went,Total)
>
> glm<-glm(Went~fTemp+fSal+fTemp*fSal,family="binomial",data=DF)
>
> require(graphics)
>
> plot(c(5,15),c(0,1),type="n",xlab="Temperature",ylab="Probability of
> going")
>
> text(Temp,Go/Total,as.character(Sal))
>
> ld<-(seq(5,15,1))
>
>
> lines(ld,predict(glm,data.frame(fTemp=ld,fSal=factor(rep("34",length(ld)),levels=levels(fSal))),type="response"))
>
>
>
> Thank you very much for your time and expertise!
>
> Kindly,
>
> Chad
>
> [[alternative HTML version deleted]]
>
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
