R help - Oct 2010 - Plot for Binomial GLM

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






-- 
View this message in context:
http://r.789695.n4.nabble.com/Plot-for-Binomial-GLM-tp2954406p2954406.html
Sent from the R help mailing list archive at Nabble.com.

On Mon, Oct 4, 2010 at 7:21 AM, klsk89 <karenklsk89 at yahoo.com>
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.:)
What about xtabs?   For instance:

xtabs(deadalive ~ Dose + Sex, data = rat.toxic)

Regarding graphs, take a look at faceting in ggplot2 (or lattice).
You can get something close to the 3 way table but in graphical form
that way.  I am not sure if this is completely up and running yet, but
I know there has been work linking ggobi with R.  I have seen a few
demonstrations that looked quite promising, and it may work well for
you to visualize three variables at once (and interactively).  Here is
the link:  http://www.ggobi.org/rggobi/
>
> Here's my data.
>
>>
rat.toxic<-read.table(file="Rats.csv",header=T,row.names=NULL,sep=",")
>> attach(rat.toxic)
      ^  why attach it?
>> 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
Please tell me that after this, you converted the counts of dead and
alive into a single variable that had a 0 or 1 if dead and the
opposite as alive before you used it as the dependent variable in your
logistic regression.
> 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
>
>
>
>
>
>
> --
> View this message in context:
http://r.789695.n4.nabble.com/Plot-for-Binomial-GLM-tp2954406p2954406.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at 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.
>


-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/

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