[R] Plot for Binomial GLM
Joshua Wiley
jwiley.psych at gmail.com
Mon Oct 4 21:08:00 CEST 2010
Hi,
Dennis was kind of enough to remind me that glm() can take a two
column matrix, which is probably what you did with deadalive. He also
gave a rather elegant graphing solution using xyplot:
xyplot(Alive/20 ~ Dose, data = rat.toxic, groups = Sex, type = c('p', 'a'))
Josh
On Mon, Oct 4, 2010 at 8:23 AM, Joshua Wiley <jwiley.psych at gmail.com> wrote:
> 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/
More information about the R-help
mailing list