# 6. Harjoitukset

# H6.1
ldose <- rep(0:5, 2)
numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)
number <-rep(20,12) # painomuuttuja
y <- numdead/number # vastemuuttuja
sex <- factor(rep(c("M", "F"), c(6, 6)))
SF <- cbind(numdead, numalive=20-numdead)

budworm.lg <- glm(SF ~ sex*ldose,weight=number, family=binomial)

summary(budworm.lg, cor=F)

#Coefficients:
#            Estimate Std. Error z value Pr(>|z|)    
#(Intercept)  -2.9935     0.5525  -5.418 6.02e-08 ***
#sexM          0.1750     0.7781   0.225    0.822    
#ldose         0.9060     0.1671   5.424 5.84e-08 ***
#sexM.ldose    0.3529     0.2699   1.307    0.191    
#---
lkdose<-ldose-3
 lkdose
# [1] -3 -2 -1  0  1  2 -3 -2 -1  0  1  2
budworm.lg3 <- glm(SF ~ sex*lkdose, family=binomial)

budworm.lg3 <- update(budworm.lg, . ~ sex*I(ldose-3))
summary(budworm.lg3, cor=F)$coefficients

#                    Estimate Std. Error   z value     Pr(>|z|)
#(Intercept)       -0.2754324  0.2304854 -1.195010 2.320831e-01
#sexM               1.2337256  0.3769385  3.273016 1.064064e-03
#I(ldose - 3)       0.9060363  0.1670511  5.423707 5.837539e-08
#sexM.I(ldose - 3)  0.3529131  0.2699404  1.307374 1.910856e-01

# 6.2
##########################
Onnettomuudet Cambridgessa 1978-81
Altham, s. 39
##########################

# y = onnettomuuksien lkm
y <- scan()
11 9 4 4 20 4
 
# v = arvioitu liikennemäärä
v <- scan()
2206 3276 1999 1399 2276 1417
 
rd <- c(1,1,1,2,2,2) 
# rd <- c("TR","TR","TR","MR","MR","MR")
# rd = 1 on Trumpington rd, rd = 2 on Mill rd
ToD <- c(1,2,3,1,2,3); ToD <- rep(c(1,2,3),2)
# ToD = 1,2,3 tarkoittaa aikoja 7-9.30, 9.30-15, 15-18.30.
# Huom. että muuttujan nimen "t" käyttö saattaa johtaa viheilmoitukseen,
# koska t() on R:n funktio (transpoosi)

RD <- factor(rd); ToD <- factor(ToD)

#####
# (a)
#####
acc <- glm(y ~ RD + ToD, family = poisson)
acc$fitted.values
        1         2         3         4         5         6 
# 6.923077 13.384615  3.692308  8.076923 15.615385  4.307692 
acc$residuals
          1           2           3           4           5           6 
# 0.58885263 -0.32758689  0.08333355 -0.50472818  0.28078730 -0.07142838 
y
#[1] 11  9  4  4 20  4
y-acc$fitted.values
         1          2          3          4          5          6 
# 4.0769231 -4.3846154  0.3076923 -4.0769231  4.3846154 -0.3076923 

residuals(acc,type="response")
residuals(acc,type="deviance")

acc1 <- glm(y ~ RD*ToD, family = poisson)

# Kaksi aikatyyppiä: ruuhka ja tavallinen
# Tämä
##### Rush, 2x3-taulukko ja 2x2-taulukko: 2 tyyliä
# 1. 
rush <- c(1,2,1,1,2,1); RUSH <- factor(rush)
acc.rush<- glm(y ~ RD + RUSH, family = poisson)
acc.rush$fitted.values

acc.rush$residuals

y-acc.rush$fitted.values


# 2. tyyli 2x2-taulukko:
# 2.
y <- scan()
15 9 8 20
 
rd <- c(1,1,2,2); RD<-factor(rd)
rush1 <- c(1,2,1,2); RUSH1 <- factor(rush1)
acc.rush2 <- glm(y ~ RD + RUSH1, family = poisson)
acc.rush2
##############

# 6.2 (b) 
# Tie, vuorokauden aika ja liikennemäärät
##############
lv <- log(v); 
accidents <- glm(y ~ RD + ToD + lv, family = poisson)
accidents

#summary(accidents)
rush <- c(1,2,1,1,2,1); RUSH <- factor(rush)
accidents.rush <- glm(y ~ RD + RUSH + lv, family = poisson)
accidents$deviance
accidents.rush$deviance
#[1] 1.882272
#[1] 1.889615

pchisq(1.882272,df=1, ncp=0, lower.tail = F)
pchisq(1.889615,df=2, ncp=0, lower.tail = F)
# 0.1700761
# [1] 0.3887544
anova(accidents,test="Chisq")	# Analysis of Deviance Table
anova(accidents.rush,test="Chisq")
# anova antaa deviannsianalyysin taulukon

# 6.3 (b) SUE:n kovarianssimatriisi

y <- scan()
0  0  3  0  1  1  1  2  2
4  2  8  0  3  4  5  2  2
2  5  4  3 15 12  7 14 6
10 14  8 19 10  7 20 10 19
 
i <- 1:36
aids.reg <- glm(y ~ i,family=poisson)

# Kovarianssimatriisi  
summary(aids.reg, cor=T)
# s12= -0.9492*0.211933*0.007708
#[1] -0.001550594
 0.211933^2
#[1] 0.0449156
0.007708^2
# [1] 5.941326e-05

###########

aids.reg$coefficients 
# 0.03966095  0.07956924 
e<-exp(-0.03966095 + 0.07956924*i)
a<-sum(e); b<-sum(i*e); c<-sum(i^2*e)
Im<-array(c(a,b,b,c),dim=c(2,2))
Im<- matrix(c(a,b,b,c),2)

#Im
#          [,1]       [,2]
#[1,]  207.8421   5425.139
#[2,] 5425.1393 157148.744
#solve(Im)
             [,1]          [,2]
#[1,]  0.048653134 -1.679619e-03
#[2,] -0.001679619  6.434775e-05


# 6.7
###############
# Royal Ascot ja Henley Regatta
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
#######
# (a)
saturated <- glm(r ~ ROW*COL, family=poisson)
independence <- glm(r ~ ROW+COL, family=poisson)
saturated;independence
independence$deviance
#[1] 0.677331
pchisq(0.677331,df=1, ncp=0, lower.tail = F)
# [1] 0.4105073 ; Malli sopii hyvin

#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) Yhdysvaikutustermi ei ole merkitsevä, joten riippumattomuusmalli jää voimaan
# (c) Malli "independent" hyväksytään, kiinnijäämisriskissä ei eroa.
#############