The R code for this model and other models in this section for the boy scout example is in the R program scout.R. Here is part of the code that works with the \(2\times2\) (marginal) table between scouting and delinquency.

#### corresponds to the data scout1 in scout.SAS
#### Fits a Logistic regression with S="scout" or "nonscout"

S=factor(c("scout","nonscout"))
Sscout=(S=="scout")
Snonscout=(S=="nonscout")
y=c(33,64)
n=c(376,424)
count=cbind(y,n-y)
result=glm(count~Sscout+Snonscout,family=binomial("logit"))
summary(result)

Part of the R output includes:

Coefficients: (1 not defined because of singularities)
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -1.7272     0.1357 -12.732  < 2e-16 ***
SscoutTRUE     -0.6140     0.2272  -2.702  0.00688 ** 
SnonscoutTRUE       NA         NA      NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Note that R sets the baseline to "non-scout" because it comes before "scout" in alphabetical order. The estimated coefficient for scout is

\(\hat{\beta}_1=−0.6140\)

and is identical to the log-odds ratio from the analysis of the \(2\times2\) table. The standard error for \(\hat{\beta}_1\), 0.2272, is identical to the standard error that came from the \(2\times2\) table analysis.

Also, in this model, \(\beta_0\) is the log-odds of delinquency for nonscouts (\(x_i=0\)). Looking at the \(2\times2\) table, the estimated log-odds for non-scouts is

\(\log(64/360) = −1.7272\)

which agrees with \(\hat{\beta}_0\) from the logistic model.

The goodness-of-fit statistics are shown below:

   Null deviance: 7.6126e+00  on 1  degrees of freedom
Residual deviance: 4.0190e-14  on 0  degrees of freedom
AIC: 15.083

As we have seen before, the residual deviance is essentially 0 because the model is saturated. The null deviance corresponds to the deviance goodness-of-fit statistics in the SAS output. Here, the deviance \(G^2 = 7.6126\) is identical to the \(G^2\) for testing independence in the \(2\times2\) table.

Here is how we can fit this model in R (from scout.R, corresponds to the scout2 data in the scout.sas program).

S=factor(c("low","medium","high"))
y=c(53,34,10)
n=c(265,270,265)
count=cbind(y,n-y)
Smedium=(S=="medium")
Shigh=(S=="high")
result=glm(count~Smedium+Shigh,family=binomial("logit"))
summary(result)

Notice that we only use “Smedium” and “Shigh” in the model statement in glm(). So we set the baseline as "low" for this model.

R output:

Call:
glm(formula = count ~ Smedium + Shigh, family = binomial("logit"))

Deviance Residuals: 
[1]  0  0  0

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.3863     0.1536  -9.027  < 2e-16 ***
SmediumTRUE  -0.5512     0.2392  -2.304   0.0212 *  
ShighTRUE    -1.8524     0.3571  -5.188 2.13e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 3.6252e+01  on 2  degrees of freedom
Residual deviance: 6.8390e-14  on 0  degrees of freedom
AIC: 20.942

Number of Fisher Scoring iterations: 3

The null deviance is the \(G^2\) that corresponds to the deviance goodness-of-fit statistics found in the SAS output. Here, the deviance \(G^2 = 36.252\), so we can conclude that the delinquency is related to socioeconomic status. The test of the global null hypothesis \(\beta_1=\beta_2=0\) is equivalent to the usual test for independence in the \(3\times2\) table. The estimated coefficients and SEs are as we predicted, and the estimated odds ratios are

\(\exp(−.5512) = 0.576 = 1/1.735\),
\(\exp(−1.852) = 0.157 = 1/6.375\).