# TYÖARKKI 10. # # # Harjoitellaan logaritmisesti lineaarista mallintamista # Poissonin jakauman avulla. # sink()-komennon käyttö tulostuksen ohjaamisessa. # Selvitä komentojen tarkoitus ja tulkitse tulostus. # # Aineistona uusien rekisteröityjen AIDS-tapausten määrä # marraskuuhun 1985 mennessä (A. Sykes 1986). Mallinnetaan # tapausten lukumäärä kuukauden järjestysnumeron avulla. # 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 plot(i,y) #curve aids.reg <- glm(y ~ i,family=poisson) # poisson pienellä vaikka nimi plot(i,y) f <- aids.reg$fitted.values lines(f) summary(aids.reg) Call: glm(formula = y ~ i, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -2.4196 -1.1553 -0.2742 0.7264 2.8500 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.039661 0.211933 0.187 0.852 i 0.079569 0.007708 10.323 <2e-16 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 190.17 on 35 degrees of freedom Residual deviance: 62.36 on 34 degrees of freedom AIC: 177.69 Number of Fisher Scoring iterations: 4 ################ dchisq(x, df, ncp=0, log = FALSE) pchisq(q, df, ncp=0, lower.tail = TRUE, log.p = FALSE) qchisq(p, df, ncp=0, lower.tail = TRUE, log.p = FALSE) rchisq(n, df) pchisq(62.36, 34, ncp=0, lower.tail = F) #[1] 0.002131935 ################# mu <- exp(0.039661+0.079569*i) plot(i,mu) ################# #Lasketaan devianssi suoraan kaavalla #1.7(b) ye <- y/aids.reg$fitted.values ye1 <- ye[ye>0] y1 <- y[y>0] 2*sum(y1*log(ye1)) # [1] 62.35956 ################# #1.7(c) yy <- scan() 5 5 6 8 7 5 6 5 7 15 12 7 14 6 10 14 8 19 10 7 20 10 19 ii <- c(3.5,8,10.5,12,14,16,18,20,21.5,23:36) aids1.reg <- glm(yy ~ ii,family=poisson) # poisson pienellä vaikka nimi summary(aids.reg) Call: glm(formula = yy ~ ii, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -1.9116 -0.8964 -0.0917 0.7646 1.7659 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.349121 0.226333 5.961 2.51e-09 *** ii 0.037631 0.008294 4.537 5.70e-06 *** --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 48.549 on 22 degrees of freedom Residual deviance: 26.155 on 21 degrees of freedom AIC: 122.92 ############### pchisq(26.155, 22, ncp=0, lower.tail = F) #[1] 0.2450866 ############### #TEHTÄVÄ 1.8 # Lasketaan devianssi us:n kautta ######### # uf:n maksimi w.f:n vallitessa Lmaxf <- prod(dpois(y, lambda=y)) -2*sum(log(dpois(y, lambda=y))) # uf:n maksimi w.c:n vallitessa mu.hat <- aids.reg$fitted.values Lmaxc <- prod(dpois(y, lambda=mu.hat)) #Lmaxc #[1] 1.920280e-38 # uf:n maksimien suhde Rn <- Lmaxf/Lmaxc # Rn # 3.477066e+13 2*log(Rn) #[1] 62.35959 Lf <- log(Lmaxf) Lc <- log(Lmaxc) c(Lf,Lc,2*Lf,2*Lc,2*Lf-2*Lc) #################### Lf <- sum(log(dpois(y, lambda=y))) Lc<- sum(log(dpois(y, lambda=mu.hat))) -2*Lc; -2*Lf; -2*Lc-(-2*Lf) #################### # # # TEHTÄVÄ 1.12 # #################### aids.reg$weights aids.reg$fitted.values summary(aids.reg, correlation = T) # 0.211933*0.007708*-0.9492 # kovarianssi = -0.001550594 w <- aids.reg$weights wx <- w*i X <- matrix(c(rep(1,36),1:36),ncol=2,byrow=F) wwx <- rbind(w,wx) XWX <- wwx%*%X solve(XWX) # [,1] [,2] #[1,] 0.044915453 -0.0015506634 #[2,] -0.001550663 0.0000594131 sqrt(0.044915453);sqrt(0.0000594131) # 0.2119327 0.00770799 -0.0015506634/(sqrt(0.044915453)*sqrt(0.0000594131)) # -0.9492455 #################### #sink("C:\\Kurssit\\TMallit\\Rscripts\\temp.txt") # Tämän jälkeen tulostus ohjautuu tiedostoon # C:\\Kurssit\\TMallit\\Rscripts\\temp.txt # sink() palauttaa tulostuksen näytölle. aids.reg # Ei tulostusta näytölle summary(aids.reg) # Ei tulostusta näytölle q() # Edellisessä esimerkissä devianssi on suuri verrattuna Khi2_34 # muuttujan odotusarvoon, mutta aineistossa on monia pieniä # frekvenssejä (< 5). Khi2-approksimaatio saattaa olla siis huono. # Tilannetta voidaan parantaa yhdistämällä viereisiä soluja. # Tee tällainen uusi analyysi ja vertaile tuloksia # # 2. aineisto # The Independent julkaisi 10.11.1999 artikkelin # "GM trees pose risk of ´superweed´ calamity" # Artikkelissa esitetään aineisto, jossa kerrotaan kenttäkokeiden # yhteydessä luontoon päässeiden GM puulajikkeiden lukumäärä # vuosona 1988 - 1998. # # 1988 89 90 91 92 93 94 95 96 97 98 # 1 1 0 2 2 2 1 1 6 3 5 # # Esimerkiksi 1988 oli yksi GM lajike (eurooppalainen haapa) # oli päässyt kenttäkokeista luontoon, mutta 1998 5 lajiketta. # Kaikkiaan 24 lajiketta oli päässyt luontoon. Aineisto on osa # 17 maassa suoritetuista ainakin 116 GM puilla tehdystä kokeesta, # joissa mukana 24 lajiketta. # Kokeile aineistoon Poissonin regressiota. # # 3. aineisto # The Independent julkaisi 18. lokakuuta 1995 aineiston, jossa # oli annettu ministerien ministerien erot tehtävistään kesken virkakauden # seuraavista syistä: seksiskandaali, talousskandaali, epäonnistuminen, # poliittinen periaate tai julkinen arvostelu. Aineisto alkaa työväenpuolueen # hallituksella 1945-51. # # 45-51 51-55 55-57 57-63 63-64 64-70 70-74 74-76 76-79 79-90 90-95 # lab con con con con lab con lab lab con con # 7 1 2 7 1 5 6 5 4 14 11 # # Sex 0 0 0 2 1 0 2 0 0 1 4 # Fin 1 0 0 0 0 0 2 0 0 0 3 # Fai 2 1 0 0 0 0 0 0 0 3 0 # Pol 3 0 2 4 0 5 2 5 4 7 3 # Pub 1 0 0 1 0 0 0 0 0 3 1 # # Onko Työväenpuolueen ja Konservatiivien eroamisaste sama? Offset-muuttuja # log(t) otettava malliin. t on hallituksen aika vallassa. Ohessa on eritelty # erot myös syyn mukaan. ####################### ####################### # 4. Tieliikenneonnettomuudet Cambridgessa 1978-81 # Kuinka onnettomuudet riippuvat liikenteen määrästä, # tiestä ja vuorokaudenajasta? (ks. Altham s. 39) ################# # TEHTÄVÄ 1.15 # ################# # 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 = 1 on Trumpington rd, rd = 2 on Mill rd ToD <- c(1,2,3,1,2,3) # 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) lv <- log(v); RD <- factor(rd); ToD <- factor(ToD) accidents <- glm(y ~ RD + ToD + lv, family = poisson) accidents summary(accidents) pchisq(1.8823,df=2, ncp=0, lower.tail = F) # 0.3901789 anova(accidents) # Analysis of Deviance Table # anova antaa deviannsianalyysin taulukon drop.road <- update(accidents,.~. -RD) pchisq(7.591,df=2, ncp=0, lower.tail = F) # 0.02247167 drop.ToD <- update(accidents,.~. -ToD) drop.lv <- update(accidents,.~. -lv) summary(drop.road) rush <- c(1,2,1,1,2,1); RUSH <- factor(rush) acc.new <- glm(y ~ RD + RUSH + lv,poisson) # Tulkitse faktori "RUSH" summary(acc.new)