On 2010-10-04 8:21, klsk89 wrote:
Hi i would like to use some graphs or tables to explore the data and make
some sensible guesses of what to expect to see in a glm model to assess if
toxin concentration and sex have a relationship with the kill rate of rats.
But i cant seem to work it out as i have two predictor
variables~help?Thanks.:)
Here's my data.
rat.toxic<-read.table(file="Rats.csv",header=T,row.names=NULL,sep=",")
attach(rat.toxic)
names(rat.toxic)
[1] "Dose" "Sex" "Dead" "Alive"
rat.toxic
Dose Sex Dead Alive
1 10 F 1 19
2 10 M 0 20
3 20 F 4 16
4 20 M 4 16
5 30 F 9 11
6 30 M 8 12
7 40 F 13 7
8 40 M 13 7
9 50 F 18 2
10 50 M 17 3
11 60 F 20 0
12 60 M 16 4
13 10 F 3 17
14 10 M 1 19
15 20 F 2 18
16 20 M 2 18
17 30 F 10 10
18 30 M 8 12
19 40 F 14 6
20 40 M 12 8
21 50 F 16 4
22 50 M 13 7
23 60 F 18 2
24 60 M 16 4
glm2<-glm(deadalive~Dose*Sex,family=binomial,data=rat.toxic)
anova(glm2,test="Chi")
Analysis of Deviance Table
Model: binomial, link: logit
Response: deadalive
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 23 225.455
Dose 1 202.366 22 23.090<2e-16 ***
Sex 1 4.328 21 18.762 0.0375 *
Dose:Sex 1 1.149 20 17.613 0.2838
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
summary(glm2)
Call:
glm(formula = deadalive ~ Dose * Sex, family = binomial, data = rat.toxic)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.82241 -0.85632 0.06675 0.61981 1.47874
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.47939 0.46167 -7.537 4.83e-14 ***
Dose 0.10597 0.01286 8.243< 2e-16 ***
SexM 0.15501 0.63974 0.242 0.809
Dose:SexM -0.01821 0.01707 -1.067 0.286
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 225.455 on 23 degrees of freedom
Residual deviance: 17.613 on 20 degrees of freedom
AIC: 91.115
Number of Fisher Scoring iterations: 4
This is pretty much the budworm example in MASS (the book).
I would produce a plot of Prob(dead) vs dose with separate
lines for M/F:
dose <- seq(10, 60, length=51)
ypF <- predict(glm2, data.frame(Dose=dose, Sex="F"), type="response")
ypM <- predict(glm2, data.frame(Dose=dose, Sex="M"), type="response")
plot(c(10,60), c(0,1), type="n")
lines(dose, ypF, col=2)
lines(dose, ypM, col=4)
text(locator(1), "F", col=2)
text(locator(1), "M", col=4)
See recent posts by Ben Bolker for confidence bands.
-Peter Ehlers
______________________________________________
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.
