# TYÖARKKI 6.
#
# Logistinen regressio binomijakaumalle.
#   Data on kirjasta "Modelling Binary Data" (D. Collett, 1991).
# Tutkitaan lentokoneiden rakennuksessa käytettävien metallin kiinnittimien
# puristuskestävyyttä. Painekuormitus kasvaa 2500:sta 4300:aan 200:n yksikön
# (psi) välein. Aineistossa
#	n = testattavien kiinnittimien lkm kussakin kuormituksessa
#	r = särkyvien kiinnittimien lkm kussakin kuormituksessa.
# Oletetaan, että r_i ~ Binomial(n_i,pi_i), i = 1, ..., 10 ja r_1, ...,r_10 ovat
# riippumattomia. Mallinnetaan parametrien pi_i riippuvuus arvoista Load_i.
# Käytetään soveltuvin osin graafisia tarkasteluja.
# Malli on
#	log(pi_i/(1-pi_i)) = a + b x load_i.
# Tämä on glm:n logit-linkki.
# Huomaa, että riippuvana muuttujana on p = r/n ja "painoina" n.
# glm:n syntaksiin voi tutustua komennolla
# help(glm)
# Data on tiedostossa "Alloyf".
data <- read.table(("C:\\Kurssit\\TMallit\\Datat\\Alloyf.txt"),header=T)
attach(data)
data
   Load   n  r
1    25  50 10
2    27  70 17
3    29 100 30
4    31  60 21
5    33  40 18
6    35  85 43
7    37  90 54
8    39  50 33
9    41  80 60
10   43  65 51
# Huomaa, että "attach(data)" ei pyyhi yli jo olemassaolevia muuttujia.
# Jos sinulla jo on mja r, se on poistettava: rm(r)
p <- r/n
p
plot(Load,p)
ex6 <- glm(p ~ Load,weights=n,family=binomial)
names(data); summary(data)
plot(Load,p,type="l")
ex6; summary(ex6)
names(ex6)
plot(ex6,ask=T)  # Pyydetään diagnostisia kuvioita
#Hit <Returnto see next plot: 
#Hit <Returnto see next plot: 
#
# Seuraavassa vaihdellaan linkkifunktiota. Edellä on käytetty oletus-linkkiä,
# joka binomijakauman tapauksessa on logit.
ex6.1 <- glm(p ~ Load,family=binomial(link=logit),weights=n)
ex6.p <- glm(p ~ Load,family=binomial(link=probit),weights=n)
ex6.cll <- glm(p ~ Load,family=binomial(link=cloglog),weights=n)
summary(ex6.1)
summary(ex6.p)
summary(ex6.cll)  # Komplementaarinen loglog-linkki
q()
#
# AIC on Akaiken informaatiokriteeri, joka määritellään
#	AIC = -2 x logaritmoidun usk.funktion maksimi + 2 x sov. mallin 
#	      parametrien lkm
# Logaritmoidussa usk.funktiossa on mukana vain datasta riippuva osa, ja
# samoin myös AIC:ssä
# Kaikkien sovitettavien mallien param. lkm on 2. Malleja vertailtaessa
# haetaan pienintä AIC:n arvoa.
# Laske parametrien suhde a/b eri malleille. Mitä huomaat?
#   Mikä linkkifunktio antaa parhaan sovitteen, ts. pienimmän devianssin?
# Yleensä logit ja probit sopivat aineistoon jokseenkin yhtä hyvin.
# Varoitus: load on R:ssä funktio, joten on syytä käyttää mjan nimenä
# "Load" nimen "load" sijasta.