## Teht 1
laudatur <- read.table("http://users.jyu.fi/~junyblom/laudatur.dat",
header=T)
attach(laudatur)
koulu<-factor(koulu)
sp<-factor(sp)
malli4<-glm(laud/lkm~koulu+sp+I(ruotsi-8),family=binomial("logit"),weights=lkm)
summary(malli4)
kertoimet<-exp(coef(malli4))
kertoimet
kertoimet[1]
kertoimet[1]/(1+kertoimet[1])
# On koulun yksi naisten, joiden peruskoulun arvosana oli 8, laudaturin
# todennäköisyys

# Tn:t eri koulujen naisille
naiset<-kertoimet[1]*c("koulu1"=1,kertoimet[2:4])
(naiset)/(1+naiset)

# Tn:t eri koulujen miehille
miehet<-naiset*kertoimet[5]
miehet/(1+miehet)

# Oddsien suhde, kun numero eroaa yhdellä
kertoimet[6]

summary(malli4)
1-pchisq(15.612,df=33)
abline(v=15.6)
# On riittävä

malli40<-glm(laud/lkm~sp+I(ruotsi-8),family=binomial("logit"),weights=lkm)
anova(malli40,malli4,test="Chisq")
# Koulut poikkeavat merkitsevästi toisistaan



## Teht 2

plot(malli4$lin, log(laud/(lkm-laud)), xlab="Lineaarinen prediktori",
ylab="Empiirinen logit")
abline(0,1)

plot(malli4$lin, laud/lkm, xlab="Lineaarinen prediktori",
ylab="Todennäköisyys")
lines(sort(malli4$lin), fitted(malli4)[order(malli4$lin)])


## Teht 3

malli4$fitted*lkm
# Monet ovat alle 1:n

devianssit<-vector(length=1000)

for(i in 1:1000) {
sim<-rbinom(39,lkm,malli4$fitted)
devianssit[i]<-glm(sim/lkm~koulu+sp+I(ruotsi-8),family=binomial("logit"),
weights=lkm)$deviance
}

hist(devianssit,breaks=25, prob=TRUE)
curve(dchisq(x, df=33), from=0, to=60, add=TRUE)
# odotetut freqvenssit
lkm*fitted(malli4)
mean(devianssit)
var(devianssit)
# Ei ehkä noudata Khi^2(33) - jakaumaa

# Tehtävä 4

licence <- read.table("http://users.jyu.fi/~junyblom/licence.dat", header=TRUE)
attach(licence)

sex <- factor(sex)
age <- factor(age)
training <- factor(training)

licence.glm <- glm(fail/total ~ sex + age + training, weights=total, family=binomial)
summary(licence.glm)

1-pchisq(3.0798, 7)
# malli riittävä

# vakio: tn, että opetusluvalla harjoitellut 18-19 vuotias nainen epäonnistuu ajokokeessa
# vakio: exp(0.19055)/(1+exp(0.19055)) = 0.5474939

round(exp(coef(licence.glm)[-1]),2)

# miehillä 0.55-kertainen riski verrattuna naisiin
# autokoulun käyneillä 0.45-kertainen riski verrattuna opetusluvalla harjoitelleisiin

# kun verrataan samaa sukupuolta ja samalla menetelmällä ajamaan oppineita
# 20-29 vuotiailla on 1.23-kertainen riski epäonnistua ajokokeessa verrattuna 18-19 vuotiaisiin
# 30-39 vuotiailla 0.93-kertainen riski
# yli 40 vuotiailla 3.76-kertainen riski

t(round(exp(confint(licence.glm)[-1,]),2))

# yli 40 vuotiailla huomattavasti suurempi riski epäonnistua ajokokeessa

licence.glm0 <- glm(fail/total ~ sex + training, weights=total, family=binomial)
anova(licence.glm0, licence.glm, test="Chisq")
# ikäryhmät poikkeavat toisistaan merkitsevästi

# Tehtävä 5

# 18-19 vuotiaat autokoulun käyneet naiset 
# 0.19055-0.79113=-0.60058
tn1a <- exp(-0.60058)/(1+exp(-0.60058)) 
tn1a

# 20-29 vuotiaat autokoulun käyneet naiset
# 0.19055-0.79113+0.20939=-0.39119
tn2a <- exp(-0.39119)/(1+exp(-0.39119))
tn2a

# 30-39 vuotiaat autokoulun käyneet naiset
# 0.19055-0.79113-0.07373=-0.67431
tn3a <- exp(-0.67431)/(1+exp(-0.67431))
tn3a

# yli 40 vuotiaat autokoulun käyneet naiset
# 0.19055-0.79113+1.32365=0.72307
tn4a <- exp(0.72307)/(1+exp(0.72307))
tn4a

# 18-19 vuotiaat autokoulun käyneet miehet 
# 0.19055-0.79113-0.59649=-1.19707
tn1b <- exp(-1.19707)/(1+exp(-1.19707))
tn1b

# 20-29 vuotiaat autokoulun käyneet miehet
# 0.19055-0.79113-0.59649
tn2b <- exp(-0.98768)/(1+exp(-0.98768))
tn2b

# 30-39 vuotiaat autokoulun käyneet miehet
tn3b <- exp(-1.2708)/(1+exp(-1.2708))
tn3b

# yli 40 vuotiaat autokoulun käyneet miehet
tn4b <- exp(0.12658)/(1+exp(0.12658))
tn4b

tn1a/tn1b
tn2a/tn2b
tn3a/tn3b
tn4a/tn4b


######## Tehtävä 7 ########

framingham <- read.table("http://users.jyu.fi/~junyblom/framingham.dat", header=T)
osadata <- framingham[framingham$period=="0-9",]
attach(osadata)

sex <- factor(sex)
sbp <- factor(sbp)
age <- factor(age)

sairaus.malli <- glm(chd ~ age+sex+sbp, family=poisson, offset=log(1e-3*pyears))
summary(sairaus.malli)

### (a) ###

exp(1.0010)

# Naisten riski sairastua noin 2.72-kertainen verrattuna miehiin


### (b) ###

exp(1.1662)

# Niillä, joilla on korkea verenpaine, on noin 3.2-kertainen riski sairastua


### (c) ###

exp(0.2735)		# 50-54-vuotiaat
exp(1.2244)		# 55-59-vuotiaat	
exp(1.3519)		# 60-62-vuotiaat

detach(framingham)
detach()