Statystyka, prognozowanie, ekonometria, data mining

bulva · Podporucznik **Pomógł:** 6 razy Posty: 205 Skąd: Zgierz

Witam,

Chciałbym napisać w R kod który przeprowadzi test permutacji oparty na statystyce F (z one-way ANOVA) oraz analize post-hoc opartą na HSD Tukeya w przypadku tzw block-design (czyli ma uwzgledniac kazda mozliwa kombinację bloków/warunków i sposobow leczenia )

Mam koncepcyjny problem z tym ostatnim, tj jak zaimplementować HSD Tukeya do wyników mojego testu? Jak powinny wygladac porówania? Czy tez powinny byc oparte o F statystyke czy o T? Poniżej kod jak do tej pory:

Kod:

lipids=data.frame(LLR=c( 0.73, 0.86, 0.94, 1.40, 1.62, 0.67,0.75,0.81,1.32,1.41,0.15,0.21,0.26,0.75,0.78),
Block=rep(c(1,2,3,4,5),3), Treatment=c(rep(1,5), rep(2,5), rep(3,5)) )
library(sqldf)
lipids=sqldf("select * from lipids order by Block")
lipids$Block=as.factor(lipids$Block)
lipids$Treatment=as.factor(lipids$Treatment)
lipids

silnia = function(n){if (n<= 0) 1 else n * silnia(n-1)}

N=silnia(3)^5
F_vector = rep(NA, N)
perm.data = lipids

for(i in 1:N)
{
set.seed(i)
perm.labels = c( sample(c(1,2,3),size=length(c(1,2,3)),replace=F), sample(c(1,2,3),size=length(c(1,2,3)),replace=F),
sample(c(1,2,3),size=length(c(1,2,3)),replace=F), sample(c(1,2,3),size=length(c(1,2,3)),replace=F),
sample(c(1,2,3),size=length(c(1,2,3)),replace=F) )
perm.data$Treatment=perm.labels
F_vector[i]=summary(aov(LLR~Block+Treatment,data=perm.data))[[1]][2,4]
}

F_obs=summary(aov(LLR~Block+Treatment,data=lipids))[[1]][2,4]
p_val=mean(F_vector>=F_obs)
p_val

Dane:

Kod:

# Block extreme low fairly low moderately low
# ages 15-24 0.73 0.67 0.15
# ages 25-34 0.86 0.75 0.21
# ages 35-44 0.94 0.81 0.26
# ages 45-54 1.40 1.32 0.75
# ages 55-64 1.62 1.41 0.78

Google

Crunchy · Wysłany: 2009-10-26, 22:54

Zobacz co jest w pakiecie coin.

bulva · Podporucznik **Pomógł:** 6 razy Posty: 205 Skąd: Zgierz

Witam,
Zan ten pakiet, ale tak jak napisalem jestem zainteresownay raczej napisaniem wlasnego kodu niz korzystanie z gotowego rozwiazania, chcialbym potem wykorzystac go do symulacji. Pakietem tym moge sobie wyniki sprawdzić najwyzej i porównać z moimi co też jest dobre :)

Crunchy · Wysłany: 2009-10-27, 20:45

No to musisz się pomęczyć - się chwali ;-)

Kod:

?ptukey

bulva · Podporucznik **Pomógł:** 6 razy Posty: 205 Skąd: Zgierz

Mam takie cuś:

Kod:

silnia = function(n){if (n<= 0) 1 else n * silnia(n-1)}

N=silnia(3)^5
F_vector = rep(NA, N)
Q=rep(NA,N)
perm.data = lipids
a=c(1,2,3)

for(i in 1:N)
{
set.seed(i)
perm.labels = c( sample(a, size=length(a),replace=F), sample(a,size=length(a),replace=F),
sample(a,size=length(a),replace=F), sample(a,size=length(a),replace=F),
sample(a,size=length(a),replace=F) )
perm.data$Treatment=perm.labels
F_vector[i]=summary(aov(LLR~Block+Treatment,data=perm.data, alpha=0.1))[[1]][2,4]
Q[i]=max( abs(with(perm.data, t.test(LLR[Treatment=='1'] - LLR[Treatment=='2'], data=perm.data, alpha=0.1)$statistic)),
abs(with(perm.data, t.test(LLR[Treatment=='1'] - LLR[Treatment=='3'], data=perm.data, alpha=0.1)$statistic)),
abs(with(perm.data, t.test(LLR[Treatment=='2'] - LLR[Treatment=='3'], data=perm.data, alpha=0.1)$statistic) ) )
}

T12_obs=with(lipids, t.test(lipids$LLR[Treatment=='1'] - lipids$LLR[Treatment=='2'], data=lipids, alpha=0.1)$statistic)
T13_obs=with(lipids, t.test(lipids$LLR[Treatment=='1'] - lipids$LLR[Treatment=='3'], data=lipids, alpha=0.1)$statistic)
T23_obs=with(lipids, t.test(lipids$LLR[Treatment=='2'] - lipids$LLR[Treatment=='3'], data=lipids, alpha=0.1)$statistic)

sort(Q)[N*0.90]
ifelse(T12_obs> sort(Q)[N*0.90],"Conclude the alternative", "Insuffiecient evidence to exclude null hypothesis" )
ifelse(T13_obs> sort(Q)[N*0.90],"Conclude the alternative", "Insuffiecient evidence to exclude null hypothesis")
ifelse(T23_obs> sort(Q)[N*0.90],"Conclude the alternative", "Insuffiecient evidence to exclude null hypothesis")

Wektor Q przechowuje maksimum z "parowanych" statystyk T testów studenta, potem z jego punktem krytycznym (na poziomie istotności 0.1) porównuje statystyki "nie-permutowanych" testów. Czy taka analiza post-hoc jest poprawna w przypadku eksperymentu blokowego?