[R] splitting a factor in an analysis of deviance table (negative binomial model)
Walmes Zeviani
walmeszeviani at hotmail.com
Sat May 29 01:07:29 CEST 2010
Sorry by the delay. You could do:
> my.data <- expand.grid(A=factor(1:4), B=factor(1:4), rep=1:4)
> my.data$y <- rbinom(my.data$A, 10, 0.5)
>
> model <- glm(cbind(y, 10-y)~A*B, family=binomial, data=my.data)
> anova(model, test="Chisq" )
Analysis of Deviance Table
Model: binomial, link: logit
Response: cbind(y, 10 - y)
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 63 53.528
A 3 1.6792 60 51.849 0.6416
B 3 2.9380 57 48.911 0.4013
A:B 9 13.8380 48 35.073 0.1282
>
> X <- model.matrix(~A/B, my.data)
>
> A <- X[,2:4]
> B.A1 <- X[,grep("A1:", colnames(X))]
> B.A2 <- X[,grep("A2:", colnames(X))]
> B.A3 <- X[,grep("A3:", colnames(X))]
> B.A4 <- X[,grep("A4:", colnames(X))]
>
> model2 <- glm(cbind(y, 10-y)~A+B.A1+B.A2+B.A3+B.A4, family=binomial,
> data=my.data)
> anova(model2, test="Chisq")
Analysis of Deviance Table
Model: binomial, link: logit
Response: cbind(y, 10 - y)
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 63 53.528
A 3 1.6792 60 51.849 0.64156
B.A1 3 1.2924 57 50.556 0.73093
B.A2 3 4.3126 54 46.244 0.22962
B.A3 3 1.8865 51 44.357 0.59629
B.A4 3 9.2844 48 35.073 0.02574 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
Walmes.
-----
..ooo0
...................................................................................................
..(....)... 0ooo... Walmes Zeviani
...\..(.....(.....)... Master in Statistics and Agricultural
Experimentation
....\_)..... )../.... walmeszeviani at hotmail.com, Lavras - MG, Brasil
............
(_/............................................................................................
--
View this message in context: http://r.789695.n4.nabble.com/splitting-a-factor-in-an-analysis-of-deviance-table-negative-binomial-model-tp1017744p2235345.html
Sent from the R help mailing list archive at Nabble.com.
More information about the R-help
mailing list