R-hjelp til oppgave 2 ST301 eksamen V2000

# Les dataene inn i en dataramme, gi navn til variablene i datarammen og se p? dataene:

biller<-read.table("http://www.math.uio.no/avdc/kurs/ST301/data/beetle.dat")

names(biller)<-c("dose","ant","dode")

biller

# Kontroller at dataene svarer til de som er gitt p? side 69 i BS,

# (men merk at der er det oppgitt 10-logaritmen av dosene)

# Gj?r variablene i datarammen tilgjengelige:

attach(biller)

# Finner andel d?de ved de ulike dosene og plotter dem mot log-dose

andel<-dode/ant

plot(log(dose,10), andel, pch=3, ylim=c(0,1))

# Kontroller at figuren svarer til den som er gitt p? side 69 i BS

# PUNKT a)

# Tilpass en logistisk regresjon med log-dose som kovariat og se p? resultatet.

# Siden vi her har ekte binomiske utfall, ikke bin?re, settes responsen opp

# som cbind(y,n-y) der n=antall individer i en gruppe og y=antall "suksesser" i gruppen.

biller.a<-glm(cbind(dode,ant-dode)~log(dose,10),family=binomial)

summary(biller.a)

# PUNKT b)

# En annen variant av regresjon for binomiske (og bin?re) data er ? benytte 

# en s?kalt komplement?r log-log link (se eksamensoppgaven). 

biller.b<-glm(cbind(dode,ant-dode)~log(dose,10), family=binomial(link=cloglog))

summary(biller.b)

# For ? se hvordan de to modellene passer, er det hensiktsmessig ? se p? 

# et plot av de estimerte sannsynlighetene for de to modellene og andelen d?de 

# mot log-dose

plot(log(dose,10), andel, pch=3, ylim=c(0,1))

points(log(dose,10), biller.a$fit, pch=1)

points(log(dose,10), biller.b$fit, pch=2)

legend(1.82,0.3,c("andel","modell a","modell b"),pch=c(3,1,2))

# For ? se hvilken modell som passer best, kan vi ogs? se 

# p? residualdeviansen? for de to modellene. Hva sier de deg? 

# Siden vi har binomiske data med rimelig mange biller i hver gruppe,? 

# f?r vi en formell test for modelltilpassning ved ? se sammenligne residualdeviansen

# med kji-kvadratfordelingen med 8-2=6 frihetsgrader. 

# For modellen i punkt a f?r vi P-verdien for en slik tilpassningstest ved:

1-pchisq(biller.a$dev, biller.a$df.res)

# Hvilken av de to modellene gir den beste beskrivelsen av dataene?

# Hva er problemet med den modellen som passer d?rligst?