Assignment 2: Some solutions

## Loading required package: vcd

## Loading required package: grid

## Loading required package: gnm

Exercise 2.3

data(UCBAdmissions)
UCB <- UCBAdmissions   # save typing
ftable(UCB)

##                 Dept   A   B   C   D   E   F
## Admit    Gender                             
## Admitted Male        512 353 120 138  53  22
##          Female       89  17 202 131  94  24
## Rejected Male        313 207 205 279 138 351
##          Female       19   8 391 244 299 317

1. total number of cases

sum(UCB)

## [1] 4526

2. total applicants for each dept

margin.table(UCB, 3)

## Dept
##   A   B   C   D   E   F 
## 933 585 918 792 584 714

apply(UCBAdmissions, 3, sum)

##   A   B   C   D   E   F 
## 933 585 918 792 584 714

3. proportion admitted for each dept

tab1 <- margin.table(UCB , c(1,3))
tab1

##           Dept
## Admit        A   B   C   D   E   F
##   Admitted 601 370 322 269 147  46
##   Rejected 332 215 596 523 437 668

tab1[1,] / margin.table(UCB, 3)

## Dept
##          A          B          C          D          E          F 
## 0.64415863 0.63247863 0.35076253 0.33964646 0.25171233 0.06442577

or, using prop.table

prop.table(margin.table(UCB, c(1, 3)),2)

##           Dept
## Admit               A          B          C          D          E          F
##   Admitted 0.64415863 0.63247863 0.35076253 0.33964646 0.25171233 0.06442577
##   Rejected 0.35584137 0.36752137 0.64923747 0.66035354 0.74828767 0.93557423

4. proportion in each cell admitted There are different ways to do this. Perhaps easiest is to extract the subtables of admitted vs. rejected and divide

admit <- UCB[1,,]
reject <- UCB[2,,]

admit / (admit+reject)

##         Dept
## Gender            A          B          C          D          E          F
##   Male   0.62060606 0.63035714 0.36923077 0.33093525 0.27748691 0.05898123
##   Female 0.82407407 0.68000000 0.34064081 0.34933333 0.23918575 0.07038123

Exercise 3.1: Arbuthnot data

load the data

data(Arbuthnot, package="HistData")

1. Plot Ratio vs. Year This uses old-school base R features: plot(), abline(), loess.smooth() and lines()

plot(Ratio ~ Year, 
     data=Arbuthnot, 
     type="b", 
     ylim=c(1.0,1.2), 
     ylab="M/F Ratio")
abline(h=1.0,col="red",lwd=2)
abline(h=mean(Arbuthnot$Ratio), col="blue")
Arb.sm <- loess.smooth(Arbuthnot$Year, Arbuthnot$Ratio)
lines(Arb.sm, col="blue", lwd=2)

using with() evaluates in the environment of a data frame, saves typing Arbuthnot$…

with(Arbuthnot, {
    plot(Ratio ~ Year, 
         type="b", 
         ylim=c(1.0,1.2), 
         ylab="M/F Ratio")
    abline(h=1.0, col="red", lwd=2)
    abline(h=mean(Ratio), col="blue")
    lines(loess.smooth(Year, Ratio), col="blue", lwd=2)
    })

(b) Total christenings

with(Arbuthnot, {
    Total <- Males + Females
    plot(Total ~ Year, 
         type="b", 
         ylab="Total Christenings")
    Arb.sm <- loess.smooth(Year, Total)
    lines(Arb.sm, col="blue", lwd=2)
    })

Part (a) using ggplot

library(ggplot2)
Arbuthnot |>
  ggplot(aes(x=Year, y=Ratio)) +
    geom_point() +
    geom_line() +
    geom_smooth(method = "loess", formula = y~x, color = "blue", fill="blue") +
    geom_hline(yintercept = 1, color = "red", size = 2) +
    geom_hline(aes(yintercept = mean(Ratio)), color = "lightblue", size = 1.5) +
    ylab("Male / Female Ratio")

Exercise 3.3: WomenQueue

1. Barplots

data(WomenQueue, package = "vcd")
barplot(WomenQueue, xlab="Number of Women in Queues of 10 People",
    ylab="Number of Queues", col="pink")

2. Goodness of fit

WQfit1 <- goodfit(WomenQueue, type="binomial", par=list(size=10))
unlist(WQfit1$par)

##   prob   size 
##  0.435 10.000

summary(WQfit1)

## 
##   Goodness-of-fit test for binomial distribution
## 
##                       X^2 df  P(> X^2)
## Likelihood Ratio 8.650999  8 0.3725869

3. Plot

plot(WQfit1, type="hanging")

Exercise 3.4: Saxony data

data(Saxony, package = "vcd")

1. To test goodness of fit of the binomial with an estimated value, \(\hat{p}\), only specify max. count, size as a parameter

Saxfit <- goodfit(Saxony, type="binomial", par=list(size=12))
unlist(Saxfit$par)

##      prob      size 
##  0.519215 12.000000

The summary() method gives the likelihood ratio test for goodness of fit. \(X^2 / df = 97 / 11 = 8.81\) is large (should be \(\approx 1\) for a good fitting model).

summary(Saxfit)

## 
##   Goodness-of-fit test for binomial distribution
## 
##                      X^2 df     P(> X^2)
## Likelihood Ratio 97.0065 11 6.978187e-16

Plot with a hanging rootogram

plot(Saxfit)

2. To test the model \(Bin(n=12, p = 1/2)\), you need to also specify the prob parameter

Saxfit1 <- goodfit(Saxony, type="binomial", par=list(size=12, prob=0.5))
unlist(Saxfit1$par)

## prob size 
##  0.5 12.0

Again, get the GOF test from summary() To test for additional lack of fit for the model \(Bin(n=12, p = 1/2)\) compared to the model \(Bin(n=12, p = \hat{p})\), you have to subtract manually: \(\Delta \chi^2 = 205 - 97 = 108\) on 1 df.

summary(Saxfit1)

## 
##   Goodness-of-fit test for binomial distribution
## 
##                       X^2 df     P(> X^2)
## Pearson          249.1954 12 2.013281e-46
## Likelihood Ratio 205.4060 12 2.493625e-37

Plot the constrained model with a hanging rootogram

plot(Saxfit1)