# regrese, jednoduché třídění, kontrasty, parametrizace (30. října 2008) rm(list=ls(all=TRUE)) Caj = read.csv2("data/Caj.csv") summary(Caj) attach(Caj) # plot(droga~Smes) tapply(droga,Smes,mean) tapply(droga,Smes,sd) source("errorBars.R") # variabilita v jedn. výběrech kolísá # errorBars(droga,Smes,type="SD") # úroveň a variabilita errorBars(droga,Smes,type="SE") # úroveň a PŘESNOST odhadu E Y errorBars(droga,Smes,type="CI") # úroveň a konf. interval bartlett.test(droga~Smes) # riziko heteroskedasticity # doporučená transformace n.v. s Poissonovým rozdělením plot(sqrt(droga+3/8) ~ Smes) print(prumery <- tapply(sqrt(droga+3/8),Smes,mean)) tapply(sqrt(droga+3/8),Smes,sd) errorBars(sqrt(droga+3/8),Smes,type="SD") bartlett.test(sqrt(droga+3/8)~Smes) # možný problém odstraněn # (a = lm(sqrt(droga+3/8) ~ Smes)) anova(a) anova(lm(sqrt(droga+3/8)~1)) # podmodel E Y = 1 beta_0 summary(a) # prumery prumery[-1]-prumery[1] prumery-prumery[1] # names(a) # co všechno podá a a$contrasts # informace o kontrastech # print(I <- nlevels(Smes)) # kolik úrovní faktoru Smes contr.treatment(I) # jak se z průměrů výběrů spočítají odhady coef(a): solve(cbind(1,contr.treatment(I))) # (1, C)^(-1) # # jiný základ pro contr.treatment (totožné s contr.SAS) (aTreat5 = update(a,contrasts=list(Smes=contr.treatment(n=nlevels(Smes),base=5)))) prumery mean(prumery) # tohle se ZATÍM nikde nevyskytuje anova(aTreat5) anova(a) contr.treatment(n=5,base=5) contr.SAS(I) solve(cbind(1,contr.SAS(I))) # model netřeba počítat od začátku, stačí přepočítat: update(a,contrasts=list(Smes=contr.SAS)) # (aHelmert = update(a,contrasts=list(Smes=contr.helmert))) mean(prumery) anova(aHelmert) anova(a) summary(aHelmert) contr.helmert(5) solve(cbind(1,contr.helmert(I))) # contr.sum (aSum = update(a,contrasts=list(Smes=contr.sum))) prumery mean(prumery) sum(coef(aSum)[-1]) # poslední složka odhadu a contr.sum(I) # pro uspořádané faktory (ordered) je standardně: contr.poly(I)