kahjur=read.csv("http://www.eau.ee/~ktanel/modelleerimise_koolitus_EMYs_2013/kahjur.csv", sep=";", dec=",", header=TRUE) ### Logistiline regressioon kahjur_logit1 <- glm(efekt ~ kontsentratsioon, family=binomial(logit), data=kahjur) summary(kahjur_logit1) ### Probit-regressioon kahjur_probit1 <- glm(efekt ~ kontsentratsioon, family=binomial(probit), data=kahjur) summary(kahjur_probit1) ### Kas ja kuivõrd logit- ja probit-mudelist hinnatud väärtused erinevad? x = seq(0,80) pred_logit = predict(kahjur_logit1, data.frame(kontsentratsioon=x), type="response") pred_probit = predict(kahjur_probit1, data.frame(kontsentratsioon=x), type="response") plot(x, pred_logit, type="l", xlab="Kontsentratsioon", ylab="Surevus") lines(x, pred_probit, col="red") legend(55, 0.6, c("Logit-mudel", "Probit-mudel"), col=c("black","red"), lty=1) ### Mida need mudelite parameetrid ikkagi tähendavad (ehk kuidas prognoosida)? # Logistilise regressiooni võrrand peaks eelneva analüüsi kohaselt olema kujul # logit(p) = -2.45987 + 0.07209*x, # millest suremistõenäosus 80%-se kontsentratsiooni korral on x=80 1 / (1 + exp(2.45987-0.07209*x)) # Probit-regressiooni võrrand peaks eelneva analüüsi kohaselt olema kujul # p = Fii(-1.44 + 0.414*x), kus Fii on N(0,1) jaotusfunktsioon: Fii(k)=P(Z0.5, BIN_SUGU) tabel spetsiifilisus = tabel[1,1]/sum(tabel[,1]) tundlikkus = tabel[2,2]/sum(tabel[,2]) spetsiifilisus; tundlikkus # ROC-kõver ja optimaalne tõenäosus eristamaks mehi naistest fitmees=olu_logit1$fitted.values mees_punkt=seq(min(fitmees), 1, by=0.01) nn=length(mees_punkt) tun_mees=rep(NA,nn); spe_mees=rep(NA,nn) for (i in 1:nn) { mees_progn=(fitmees>mees_punkt[i]) mees_tabel=table(mees_progn, tudeng$BIN_SUGU) tun_mees[i]=mees_tabel[2,2]/sum(mees_tabel[,2]) spe_mees[i]=mees_tabel[1,1]/sum(mees_tabel[,1]) } par(mfrow=c(1,2)) plot(1-spe_mees, tun_mees, type="b", xlab="1-Spetsiifilisus", ylab="Tundlikkus", ylim=c(0,1), xlim=c(0,1)) lines(c(0,1), c(0,1)) plot(mees_punkt, spe_mees, type="l", xlab="Piirtõenäosus", ylab="Spetsiifilisus/Tundlikkus", ylim=c(0,1), xlim=c(0,1)) lines(mees_punkt,tun_mees, lty=2) text(0.7, 0.9, "Spetsiifilisus") text(0.7, 0.4, "Tundlikkus") par(mfrow=c(1,1)) # Tundlikkus ja spetsiifilisus tõenäosuse 0,3 korral mees=predict(olu_logit1, type="response") tabel=table(mees>0.3, BIN_SUGU) tabel spetsiifilisus = tabel[1,1]/sum(tabel[,1]) tundlikkus = tabel[2,2]/sum(tabel[,2]) spetsiifilisus; tundlikkus dose.p(olu_logit1, cf = c(1,2), p = c(0.3)) # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- # Lisa. Puude raieküpsuse analüüs. puud=read.table("http://www.eau.ee/~ktanel/modelleerimise_koolitus_EMYs_2013/puud.txt", sep=";", dec=",", header=TRUE) # Moodustame binaarse tunnuse 'mature', mille väärtus on üks arenguklassi Y puhul (puu on raieküps) puud$mature = 1*(puud$ARENGUKL=="Y") attach(puud) tree_logit1 = glm(mature ~ A, family=binomial(), data=puud) summary(tree_logit1) x = seq(0,200) y = predict(tree_logit1, data.frame(A=x), type="response") plot(x, y, type="l", xlab="Vanus (aastates)", ylab="Raieküpsuse tõenäosus") points(A, mature) par(mfrow=c(2,1)) hist(A[mature==1], xlim=c(0,200)) hist(A[mature==0], xlim=c(0,200)) par(mfrow=c(1,1)) # Keerukam mudel tree_logit2 = glm(mature~A+I(A^2)+I(A^3)+I(A^4)+I(A^5),family=binomial(),data=puud) summary(tree_logit2) x = seq(0,200) y = predict(tree_logit2, data.frame(A=x), type="response") plot(x, y, type="l", xlab="Vanus (aastates)", ylab="Raieküpsuse tõenäosus") ### On's raieküpsus ja selle saavutamise vanus sõltuv ka puu liigist? tree_logit3 = glm(mature~A*factor(PE),family=binomial(),data=puud) summary(tree_logit3) library(car) Anova(tree_logit3) ### Illustreerime tulemusi kuuskede ja kaskede tarvis x = seq(0,200) yKU = predict(tree_logit3, data.frame(A=x, PE="KU"), type="response") yKS = predict(tree_logit3, data.frame(A=x, PE="KS"), type="response") plot(x, yKU, type="l", xlab="Vanus (aastates)", ylab="Raieküpsuse tõenäosus") lines(x, yKS, col="red") legend(130, 0.8, c("Kuusk", "Kask"), col=c("black","red"), lty=1) ### Raieküpsuse testi tundlikkus ja spetsiifilisus kyps=predict(tree_logit3, type="response") tabel=table(kyps>0.5, mature) tabel spetsiifilisus = tabel[1,1]/sum(tabel[,1]) tundlikkus = tabel[2,2]/sum(tabel[,2]) spetsiifilisus; tundlikkus ### Uurime edasi üksnes kuuski ja kaski kuusk = subset(puud, subset = A>0 & PE=="KU") kask = subset(puud, subset = A>0 & PE=="KS") kuusk_logit = glm(mature~A, family=binomial(), data=kuusk) kask_logit = glm(mature~A, family=binomial(), data=kask) x = seq(0,200) yKU = predict(kuusk_logit, data.frame(A=x), type="response") yKS = predict(kask_logit, data.frame(A=x), type="response") plot(x, yKU, type="l", xlab="Vanus (aastates)", ylab="Raieküpsuse tõenäosus") lines(x, yKS, col="red") legend(130, 0.8, c("Kuusk", "Kask"), col=c("black","red"), lty=1) # Raieküpsuse testi tundlikkus ja spetsiifilisus kyps_kuusk=predict(kuusk_logit, type="response") tabel_kuusk=table(kyps_kuusk>0.5, kuusk$mature); tabel_kuusk spets_kuusk = tabel_kuusk[1,1]/sum(tabel_kuusk[,1]) tund_kuusk = tabel_kuusk[2,2]/sum(tabel_kuusk[,2]) spets_kuusk; tund_kuusk kyps_kask=predict(kask_logit, type="response") tabel_kask=table(kyps_kask>0.5, kask$mature); tabel_kask spets_kask = tabel_kask[1,1]/sum(tabel_kask[,1]) tund_kask = tabel_kask[2,2]/sum(tabel_kask[,2]) spets_kask; tund_kask # ROC-kõver ja optimaalne tõenäosus eristamaks raieküpseid puid fitkuusk=kuusk_logit$fitted.values kuusk_punkt=seq(min(fitkuusk), 1, by=0.01) nn=length(kuusk_punkt) tun_kuusk=rep(NA,nn); spe_kuusk=rep(NA,nn) for (i in 1:nn) { kuusk_progn=(fitkuusk>kuusk_punkt[i]) kuusk_tabel=table(kuusk_progn, kuusk$mature) tun_kuusk[i]=kuusk_tabel[2,2]/sum(kuusk_tabel[,2]) spe_kuusk[i]=kuusk_tabel[1,1]/sum(kuusk_tabel[,1]) } par(mfrow=c(1,2)) plot(1-spe_kuusk, tun_kuusk, type="b", xlab="1-Spetsiifilisus", ylab="Tundlikkus", ylim=c(0,1), xlim=c(0,1)) lines(c(0,1), c(0,1)) plot(kuusk_punkt, spe_kuusk, type="l", xlab="Piirtõenäosus", ylab="Spetsiifilisus/Tundlikkus", ylim=c(0,1), xlim=c(0,1)) lines(kuusk_punkt,tun_kuusk, lty=2) text(0.7, 0.93, "Spetsiifilisus") text(0.75, 0.4, "Tundlikkus") par(mfrow=c(1,1)) # Tundlikkus ja spetsiifilisus tõenäosuse 0,2 korral kyps_kuusk=predict(kuusk_logit, type="response") tabel_kuusk=table(kyps_kuusk>0.2, kuusk$mature); tabel_kuusk spets_kuusk = tabel_kuusk[1,1]/sum(tabel_kuusk[,1]) tund_kuusk = tabel_kuusk[2,2]/sum(tabel_kuusk[,2]) spets_kuusk; tund_kuusk # Milline see kuuskede tarvis hinnatud logistilise regressiooni võrrand ikka oli? summary(kuusk_logit) prob=0.9 (log(prob/(1-prob))-(-19.64508))/0.22066 library(MASS) dose.p(kuusk_logit, cf = c(1,2), p = 0.9) dose.p(kuusk_logit, cf = c(1,2), p = c(0.5, 0.6, 0.7, 0.8, 0.9)) detach(puud)