# TYÖARKKI 11
#
# Käytetään Poissonin jakaumaa kaksiulotteisen kontingenssitaulukon
# logaritmisesti lineaariseen analysointiin ja verrataan tulosta vastaavaan
# binomiaaliseen formulointiin.
#   Muodostetaan 3-ulotteinen taulukko Simpsonin paradoksin
# havainnollistamiseksi
#  The Daily Telegraph (28.10.1988) julkaisi artikkelin
# "Johtajat alttiita ajamaan humalassa", jossa esitettiin 
# Royal Ascotin ja Henley Regatan yhteydessä järjestettyjen
# puhallustestien tuloksia (Royal Ascot ja Henley Regatta ovat
# tunnetusti "kosteita" englantilaisia urheilutapahtumia kuten
# hevos- ja soutukilpailut). Jääkö RA:ssa todennäköisemmin
# kiinni kuin HR:ssä (jos testataan)? RA:ssa puhallutettiin
# kaikkiaan 24+2210 ja HR:ssä 5+680. Näistä 24 positiivisia
# RA:ssa ja 5 HR:ssä.
# 
r <- scan()
24 2210   
 5 680    

# 1. rivi Royal Ascot
# 2. rivi Henley Regatta

Row <- c(1,1,2,2); Col <- c(1,2,1,2)
# Älä käytä nimeä "row", koska se on funktion nimi R:ssä
ROW <- factor(Row); COL <- factor(Col)

# Col=1 PIDÄTETTY; Col=2 EI pidätetty
##############
# Tehtävä 1.16
##############
# (a)
saturated <- glm(r ~ ROW*COL, family=poisson)
independence <- glm(r ~ ROW+COL, family=poisson)
summary(saturated)
summary(independence)
c(22.194587, 2211.805413,    6.805413,  678.194587)/2919 # solutodennäköisyydet
#0.007603490 0.757727103 0.002331419 0.232337988
# (b) Yhdysväikutustermi ei ole merkitsevä, joten riippumattomuusmalli jää voimaan
# (c) Kiinnijäämisriskissä ei eroa
#############
# khiin nelötesti 2x2-taulukoille!
#############
# Toinen tapa tarkastella asiaa
# Huomaa tietyt vastaavuudt ja selitä ne
##############
# Tehtävä 1.16
##############
# (a)
a <- c(24,2210)  # Royal Ascot
b <- c(5,680)    # Henley Regatta
tot <- a+b; p <- a/tot
Row <- c(1,2); ROW <- factor(Row)

sat <- glm(p ~ ROW,family=binomial,weights=tot)
indep <- glm(p ~ 1,family=binomial,weights=tot)
summary(sat)
summary(indep)
# sat, logit riippuu rivimuuttujasta
# indep, logit ei riipu rivimuuttujasta

# (b) Malli indep hyväksytään, kiinnijäämisriskissä ei eroa.

########################## 

# Huom. Jos X ja Y ovat riippumattomia Poissonin jakaumaa noudattavia satunnaimuuttujia,
# niin X ehdolla X+Y=n noudattaa binomijakaumaa. 

# Huom. Binomiaalisessa regressiossa ei ole välttämätöntä
# käyttää weights-määrettä, jos määritellään tehtävä hieman
# toisin.
satt <- glm(cbind(a,b) ~ ROW,binomial) 
# yhtäpitävä "sat":n kanssa
summary(satt)
indepp <- glm(cbind(a,b) ~ 1,binomial)
summary(indepp)
#########################
a1 <- c(24,5)
b1 <- c(2210,680)
satt1 <- glm(cbind(a1,b1) ~ ROW,binomial)
#
# Voi olla sangen harhaanjohtavaa summata 3-ulotteisessa
# taulukossa frekvessit jonkin dimension suhteen, esimerkiksi
# taulukossa Ascot/Henley x Pidätetty/EiPidätetty x Miehet/Naiset
# sukupuoli-dimension yli.
#
# Oletetaan esimerkiksi, että edellä esitetty 2x2-taulukko olisikin
# 3-ulotteinen taulukko (keksitty)
#  24 = 23,1  2210 = 2,2208
#   5=   3,2   680 = 340,340
# Jokaisen parin 1. luku on miesten lkm ja 2. luku naisten lkm.
#
###############
# Tehtävä 1.18
##############"
# Todetaan mallilla sat 3:n mjan yhdysvaikutus
##############
r <- scan()
23 2 1 2208
3 340 2 340
 

Row <- c(1,1,1,1,2,2,2,2); Col <- c(1,2,1,2,1,2,1,2)
sex <- c(1,1,2,2,1,1,2,2)
ROW <- factor(Row); COL <- factor(Col); SEX <- factor(sex)
sat <- glm(r ~ ROW*COL*SEX,poisson);
summary(sat)
#

#   Esimerkissä on voimakas kolmen mjan yhdysvaikutus.
# Miten pidätysaste vaihtelee Ascotin ja Henleyn välillä?
# Miten naisten pidätysaste vaihtelee Ascotin ja Henleyn välillä?
# Tämä on esimerkki Simpsonin paradoksista. Altham kutsuu esimerkin
# ilmiötä "Yulen paradoksiksi", koska G.U.Yule kirjoitti asiasta jo
# 1900-luvun alussa ja Simpson vasta n. 50 vuotta myöhemmin 1951.

# (a) Naiset

r <- scan()
1 2208   
2 340    

# 1. rivi Royal Ascot
# 2. rivi Henley Regatta

Row <- c(1,1,2,2); Col <- c(1,2,1,2)
ROW <- factor(Row); COL <- factor(Col)

# Col=1 PIDÄTETTY; Col=2 EI pidätetty

saturated <- glm(r ~ ROW*COL, family=poisson)
independence <- glm(r ~ ROW+COL, family=poisson)
summary(saturated)
summary(independence)
pchisq(4.5152, 1, ncp=0, lower.tail = F) # =0.03359496
#  Yhdysväikutustermi ei ole kovin merkitsevä, joten riippumattomuusmalli jää voimaan
#  Kiinnijäämisriskissä ei eroa
#############
# (a) Miehet

r <- scan()
23 2   
3 340    

Row <- c(1,1,2,2); Col <- c(1,2,1,2)
ROW <- factor(Row); COL <- factor(Col)

# Col=1 PIDÄTETTY; Col=2 EI pidätetty

saturated <- glm(r ~ ROW*COL, family=poisson)
independence <- glm(r ~ ROW+COL, family=poisson)
summary(saturated)
summary(independence)
pchisq(139.57, 1, ncp=0, lower.tail = F)   # = 0
#  Yhdysväikutustermi erittäin merkitsevä, joten kiinnijäämisriskeissä ei eroa