# 12. prednáška z regrese 11. listopad 2009 rm(list=ls()) # permutační výpočet p-hodnot # # připravíme všechny permutace složek Y # Y = c(-9,-11,1,19) x = c(-3,-1,1,3) (a1 = lm(Y~x)) library(e1071) (pp = permutations(4)) # matice permutací (B = nrow(pp)) # počet permutací ) dev2 = dev1 = rep(NA,B) for (i in 1:B) dev1[i] = deviance(update(a1,Y[pp[i,]]~.)) dev1 mean(dev1<=dev1[1]) # rel. četnost "nejhorších" RSS summary(a1) # a2 = lm(Y~x+I(x^2)) for (i in 1:B) dev2[i] = deviance(update(a2,Y[pp[i,]]~.)) dev2 mean(dev2<=dev2[1]) summary(a2) # Kojeni = read.csv2("data/kojeni.csv") summary(Kojeni) summary(a<- lm(hmotnost~Hoch+vek.m,data=Kojeni)) # permutací je príliš mnoho: 99! factorial(99) B = 1000; dev = rep(NA,B) # p-hodnota podmodelu E Y = konst for (i in 1:B) dev[i] = deviance(update(a,sample(hmotnost)~.)) mean(dev<=deviance(a)) # permutační p-hodnota # summary(a) psti = rep(NA,B); tStat = rep(NA,B) # permutační p-hodnota testu nulovosti koeficientu u proměnné vek.m for (i in 1:B) {summ = summary(update(a,.~.-vek.m+sample(vek.m))) psti[i] = summ$coef[3,4] tStat[i] = summ$coef[3,3] } mean(psti<=summary(a)$coef[3,4]) hist(tStat,freq=FALSE) curve(dt(x,a$df.res),-4,4,add=TRUE)# porovnání s hustotou t-rozdělení abline(v=summary(a)$coef[3,3],col="red") 1-pt(summary(a)$coef[3,3],a$df.res)# polovina p-hodnoty založené na t-rozdělení # summary(a) # permutační p-hodnota testu nulovosti koeficientu u proměnné Hoch for (i in 1:B) {summ = summary(update(a,.~.-Hoch+sample(Hoch))) psti[i] = summ$coef[3,4] tStat[i] = summ$coef[3,3] } mean(psti<=summary(a)$coef[2,4]) hist(tStat,freq=FALSE) curve(dt(x,a$df.res),-4,4,add=TRUE)# porovnání s hustotou t-rozdělení abline(v=summary(a)$coef[2,3],col="red") hist(psti) # použití erkové knihovny: library(coin) summary(a) independence_test(hmotnost~Hoch+vek.m,data=Kojeni,distribution="approximate") independence_test(hmotnost~Hoch+vek.m,data=Kojeni,distribution=approximate(B=100000))