```

References

Load packages

library(tidyverse)

Problem 1:

## Load data
data(lirat, package = "VGAM")
## A data frame with 58 observations on the following 4 variables.
## N Litter size.
## R Number of dead fetuses.
## hb Hemoglobin level.
## grp Group number. Group 1 is the untreated (low-iron) group, group 2 received injections on day 7 or day 10 only, group 3 received injections on days 0 and 7, and group 4 received injections weekly.
summary(lirat)
##        N               R                hb              grp       
##  Min.   : 1.00   Min.   : 0.000   Min.   : 3.100   Min.   :1.000  
##  1st Qu.: 9.00   1st Qu.: 0.250   1st Qu.: 4.325   1st Qu.:1.000  
##  Median :11.00   Median : 3.500   Median : 6.500   Median :1.000  
##  Mean   :10.47   Mean   : 4.603   Mean   : 7.798   Mean   :1.897  
##  3rd Qu.:13.00   3rd Qu.: 9.000   3rd Qu.:10.775   3rd Qu.:2.750  
##  Max.   :17.00   Max.   :14.000   Max.   :16.600   Max.   :4.000
## Make group 4 reference
lirat$grp <- factor(lirat$grp)
lirat$grp <- relevel(lirat$grp, ref = "4")

## Fit regular logistic model
logistic <- glm(formula = cbind(R, N - R) ~ grp + hb,
                family  = binomial(link = "logit"),
                data    = lirat)
summary(logistic)
## 
## Call:
## glm(formula = cbind(R, N - R) ~ grp + hb, family = binomial(link = "logit"), 
##     data = lirat)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.6331  -0.9684  -0.0471   1.3945   2.8256  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.1275     1.5099   0.084   0.9327  
## grp1          2.0520     1.0629   1.931   0.0535 .
## grp2         -0.4228     0.7927  -0.533   0.5938  
## grp3         -1.1006     0.9250  -1.190   0.2341  
## hb           -0.2190     0.1020  -2.148   0.0317 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 509.43  on 57  degrees of freedom
## Residual deviance: 168.91  on 53  degrees of freedom
## AIC: 250.38
## 
## Number of Fisher Scoring iterations: 5
## Add predicted values
lirat$pi_hat <- predict(logistic, type = "response")
## Standardized Pearson residuals
lirat$st_pearson_resid <- with(lirat, (R - N * p_hat) / sqrt(N * pi_hat * (1 - pi_hat)))
## Error in eval(expr, envir, enclos): object 'p_hat' not found
## Plot
ggplot(data = lirat, mapping = aes(x = N, y = st_pearson_resid)) +
    geom_point() +
    theme_bw() + theme(legend.key = element_blank())
## Error in eval(expr, envir, enclos): object 'st_pearson_resid' not found

## Overdispersion diagnostic
## Pearson X^2 / DF of residuals
sum(lirat$st_pearson_resid^2) / logistic$df.residual
## [1] 0
## Deviance / DF of residuals
logistic$deviance / logistic$df.residual
## [1] 3.186949
## Standardized Pearson residuals in rp: The vector of standardized Pearson residuals.

## Fit regular quasi-binomia model
quasibinom <- glm(formula = cbind(R, N - R) ~ grp + hb,
                  family  = quasibinomial(link = "logit"),
                  data    = lirat)
## The scale parameter is Pearson X^2 / DF of residuals
summary(quasibinom)
## 
## Call:
## glm(formula = cbind(R, N - R) ~ grp + hb, family = quasibinomial(link = "logit"), 
##     data = lirat)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.6331  -0.9684  -0.0471   1.3945   2.8256  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   0.1275     2.5889   0.049    0.961
## grp1          2.0520     1.8224   1.126    0.265
## grp2         -0.4228     1.3592  -0.311    0.757
## grp3         -1.1006     1.5860  -0.694    0.491
## hb           -0.2190     0.1749  -1.253    0.216
## 
## (Dispersion parameter for quasibinomial family taken to be 2.939931)
## 
##     Null deviance: 509.43  on 57  degrees of freedom
## Residual deviance: 168.91  on 53  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5