Reklama
|
Statystyka, prognozowanie, ekonometria, data mining
Forum miłośników statystyki - Portal Statystyczny
|
Ogłoszenie |
FORUM STATYSTYCZNE MA JUŻ 10 LAT
Znasz statystykę lub ekonometrię, metody prognozowania, data mining i chcesz pomóc w rozwoju forum statystycznego ?
Pisz na: administrator(małpa)statystycy.pl
Rozpoczął swoją działalność portal statystyczny - masz pomysł na jego rozwój ?
Drogi forumowiczu! Zanim napiszesz posta zapoznaj się z regulaminem forum i
przedstaw się
|
Przesunięty przez: mathkit 2015-03-23, 16:03 |
[R] Test permutacji |
Autor |
Wiadomość |
bulva
Podporucznik

Pomógł: 6 razy Posty: 205 Skąd: Zgierz
|
Wysłany: 2009-10-26, 17:04 [R] Test permutacji
|
|
|
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 |
|
Ostatnio zmieniony przez mathkit 2010-01-31, 10:03, w całości zmieniany 1 raz |
|
|
|
 |
Google
|
Wysłany: Reklama google.
|
|
|
|
|
|
|
Crunchy
Major Crunchy

Pomógł: 75 razy Posty: 1129 Skąd: Katowice
|
Wysłany: 2009-10-26, 22:54
|
|
|
Zobacz co jest w pakiecie coin. |
|
|
|
 |
bulva
Podporucznik

Pomógł: 6 razy Posty: 205 Skąd: Zgierz
|
Wysłany: 2009-10-27, 10:54
|
|
|
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
Major Crunchy

Pomógł: 75 razy Posty: 1129 Skąd: Katowice
|
Wysłany: 2009-10-27, 20:45
|
|
|
No to musisz się pomęczyć - się chwali
|
|
|
|
 |
bulva
Podporucznik

Pomógł: 6 razy Posty: 205 Skąd: Zgierz
|
Wysłany: 2009-10-27, 21:56
|
|
|
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? |
|
|
|
 |
|
Nie możesz pisać nowych tematów Możesz odpowiadać w tematach Nie możesz zmieniać swoich postów Nie możesz usuwać swoich postów Nie możesz głosować w ankietach Nie możesz załączać plików na tym forum Możesz ściągać załączniki na tym forum
|
Dodaj temat do Ulubionych zakładek(IE) Wersja do druku
|
| Strona wygenerowana w 0,22 sekundy. Zapytań do SQL: 19 |
|
|