> death <- read.table("death.dat",header=T)
> attach(death,pos=1)
> death
D V P count
1 w w y 53
2 w w n 414
3 w b y 0
4 w b n 16
5 b w y 11
6 b w n 37
7 b b y 4
8 b b n 139
> options(contrasts=c("contr.treatment"))
> d_v_p <- glm(count ~ D + V + P, family=poisson(link=log)) # model (D,V,P)
> summary(d_v_p)
Call: glm(formula = count ~ D + V + P, family = poisson(link = log))
Residual Deviance: 402.8353 on 4 degrees of freedom
> dvp2 <- glm(count ~ D*V + D*P + P*V, family=poisson(link=log)) # (DV,DP,PV)
> dvp2
Call:
glm(formula = count ~ D * V + D * P + P * V, family = poisson(link = log))
Degrees of Freedom: 8 Total; 1 Residual
Residual Deviance: 0.3798378
> 1 - pchisq(.3798, 1) # checking goodness of fit
[1] 0.5377103
> update(dvp2, . ~ . - D:V) # model (DP, PV)
Call:
glm(formula = count ~ D + V + P + D:P + P:V, family = poisson(link = log))
Degrees of Freedom: 8 Total; 2 Residual
Residual Deviance: 384.4273
> update(dvp2, . ~ . - D:P) # model (DV, PV)
Call:
glm(formula = count ~ D + V + P + D:V + P:V, family = poisson(link = log))
Degrees of Freedom: 8 Total; 2 Residual
Residual Deviance: 5.394039
> update(dvp2, . ~ . - P:V) # model (DV, DP)
Call:
glm(formula = count ~ D + V + P + D:V + D:P, family = poisson(link = log))
Degrees of Freedom: 8 Total; 2 Residual
Residual Deviance: 20.72975
> summary(dvp2)
Call: glm(formula = count ~ D * V + D * P + P * V, family = poisson(link = log))
Coefficients:
Value Std. Error t value
(Intercept) 4.9357837 0.08471229 58.265261
D -2.1746465 0.26375400 -8.244980
V -1.3298017 0.18478847 -7.196346
P -3.5961033 0.50662213 -7.098196
D:V 4.5949705 0.31351901 14.656114
D:P -0.8677967 0.36705684 -2.364202
P:V 2.4044427 0.60036069 4.004997
Residual Deviance: 0.3798378 on 1 degrees of freedom
Number of Fisher Scoring Iterations: 2
> pearson <- summary.lm(dvp2)$residuals # Pearson residuals
> leverage <- lm.influence(dvp2)$hat # leverage values
> adjusted <- pearson/sqrt(1 - leverage) # standardized Pearson residuals
> expected <- fitted(dvp2) # estimated expected frequencies
> cbind(count, expected, pearson, adjusted)
count expected pearson adjusted
1 53 52.8178206 0.025067201 0.4444706
2 414 414.1821794 -0.008951663 -0.4444706
3 0 0.1821795 -0.426564129 -0.4444706
4 16 15.8178207 0.045804235 0.4444706
5 11 11.1821794 -0.054477270 -0.4444706
6 37 36.8178206 0.030023748 0.4444706
7 4 3.8178233 0.093182615 0.4444706
8 139 139.1821794 -0.015442155 -0.4444706
> dv_vp <- glm(count ~ D*V + V*P, family=poisson) # log is default link
> summary(dv_vp)
Call: glm(formula = count ~ D * V + V * P, family = poisson)
Coefficients:
Value Std. Error t value
(Intercept) 4.937366 0.08458905 58.368852
D -2.190256 0.26361183 -8.308640
V -1.198864 0.16807052 -7.133100
P -3.657131 0.50632193 -7.222936
D:V 4.465384 0.30405464 14.686123
V:P 1.704547 0.52363865 3.255196
Residual Deviance: 5.394039 on 2 degrees of freedom
> pearson2 <- summary.lm(dv_vp)$residuals # Pearson residuals
> leverage2 <- lm.influence(dv_vp)$hat # leverage values
> adjusted2 <- pearson2/sqrt(1 - leverage2) # standardized Pearson residuals
> expected2 <- fitted(dv_vp) # estimated expected frequencies
> cbind(count, expected2, pearson2, adjusted2)
count expected2 pearson2 adjusted2
1 53 58.0349564 -0.66080579 -2.3122378
2 414 408.9650473 0.24897414 2.3122378
3 0 0.4025157 -0.63432196 -0.6774471
4 16 15.5974843 0.10191813 0.6774471
5 11 5.9650509 2.06037764 2.3122378
6 37 42.0349646 -0.77629580 -2.3122378
7 4 3.5974845 0.21218070 0.6774471
8 139 139.4025157 -0.03409162 -0.6774471
> def <- c(0,1,0,1); vic <- c(0,0,1,1) # now use logit models instead
> > yes <- c(53,11,0,4); no <- c(414,37,16,139)
> > binom.dat <- cbind(yes, no)
> fit.logit <- glm(binom.dat ~ def + vic, family=binomial(link=logit))
> summary(fit.logit)
Coefficients:
Value Std. Error t value
(Intercept) -2.0594573 0.1458448 -14.120879
def 0.8677967 0.3670654 2.364147
vic -2.4044429 0.6003999 -4.004735
Residual Deviance: 0.3798378 on 1 degrees of freedom
> pear.dv <- summary.lm(fit.logit)$residuals # Pearson residuals
> lev.dv <- lm.influence(fit.logit)$hat # leverage values
> adjust.dv <- pear.dv/sqrt(1 - lev.dv) # adjusted residuals
> predict.dv <- fitted(fit.logit) # estimated proportion of yes verdicts
> sample <- yes/(yes + no) # sample proportions
> cbind(sample, predict.dv, pear.dv, adjust.dv)
sample predict.dv pear.dv adjust.dv
1 0.11349036 0.11310026 0.02661769 0.4445197
2 0.22916667 0.23296207 -0.06220442 -0.4445197
3 0.00000000 0.01138622 -0.42906532 -0.4445197
4 0.02797203 0.02669806 0.09446088 0.4445197
> def <- factor(def); vic <- factor(vic)
> > cbind(def,vic)
def vic
[1,] 1 1
[2,] 2 1
[3,] 1 2
[4,] 2 2
> fit.dv <- glm(binom.dat ~ def + vic, family=binomial(link=logit))
> summary(fit.dv)
Coefficients:
Value Std. Error t value
(Intercept) -2.0594573 0.1458448 -14.120879
def 0.8677967 0.3670654 2.364147
vic -2.4044429 0.6003999 -4.004735
Residual Deviance: 0.3798378 on 1 degrees of freedom
> fit.v <- update(fit.dv, . ~ . - def) # only vic as main effect
> summary(fit.v)
Coefficients:
Value Std. Error t value
(Intercept) -1.952584 0.1335614 -14.619372
vic -1.704547 0.5236742 -3.254975
Residual Deviance: 5.394039 on 2 degrees of freedom
> pear.v <- summary.lm(fit.v)$residuals # Pearson residuals
> lev.v <- lm.influence(fit.v)$hat # leverage values
> adjust.v <- pear.v/sqrt(1 - lev.v) # adjusted residuals
> predict.v <- fitted(fit.v) # estimated proportion of yes verdicts
> cbind(sample, predict.v, pear.v, adjust.v)
sample predict.v pear.v adjust.v
1 0.11349036 0.12427185 -0.7061899 -2.3131538
2 0.22916667 0.12427185 2.2027202 2.3131538
3 0.00000000 0.02515723 -0.6425058 -0.6774973
4 0.02797203 0.02515723 0.2149161 0.6774973
>