Statystyka, prognozowanie, ekonometria, data mining Strona Główna Statystyka, prognozowanie, ekonometria, data mining
Forum miłośników statystyki - Portal Statystyczny

FAQFAQ  SzukajSzukaj  UżytkownicyUżytkownicy  GrupyGrupy  StatystykiStatystyki
RejestracjaRejestracja  ZalogujZaloguj  Chat   Regulamin  Kadra forum
PORTAL STATYSTYCZNY
 Ogłoszenie 
Rozpoczął swoją działalność portal statystyczny - masz pomysł na jego rozwój ?
Portal jest w chwili obecnej intensywnie rozwijany. Dodawane są nowe moduły, uaktualniana jest jego zawartość.
Osoby znające statystykę lub ekonometrię, metody prognozowania, data mining, a chcące pomóc w rozwoju forum statystycznego, rozbudowy portalu o dodatkowe treści, proszę o kontakt na adres e-mail: administrator(małpa)statystycy.pl

Drogi forumowiczu! Zanim napiszesz posta zapoznaj się z regulaminem forum i przedstaw się
The International Year of Statistics (Statistics2013) Free statistics help forum. Discuss statistical research, statistical consulting Smarter Poland Portal statystyczny

Poprzedni temat «» Następny temat
[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ł: 60 razy
Posty: 959
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ł: 60 razy
Posty: 959
Skąd: Katowice
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
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?
 
     
Wyświetl posty z ostatnich:   
Odpowiedz do tematu
Nie możesz pisać nowych tematów
Nie 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

Skocz do:  

Podobne Tematy
Temat Autor Forum Odpowiedzi Ostatni post
Brak nowych postów Przyklejony: Chi2 test (vs test dla dwóch wskaźników struktury)
Ania Testowanie hipotez statystycznych 103 2013-05-06, 17:21
prosie23
Brak nowych postów Przyklejony: Dokładny test Fishera (Fisher's Exact Test)
Badanie zmiennej binarnej w przekroju klas
wariancja Testowanie hipotez statystycznych 53 2013-05-13, 18:40
hamer
Brak nowych postów Przyklejony: Test istotności wahań sezonowych-test Kendalla
tomsew Testowanie hipotez statystycznych 3 2009-10-21, 18:08
Crunchy
Brak nowych postów Przyklejony: Test ADF
[Gretl]
barrel Modelowanie ekonometryczne 16 2011-05-24, 16:57
Madrich
Brak nowych postów Przyklejony: Test F
m|chał Testowanie hipotez statystycznych 32 2010-05-03, 11:00
mareczekfm

Ideą przyświecającą istnieniu forum statystycznego jest stworzenie możliwości wymiany informacji, poglądów i doświadczeń osób związanych ze statystyką, mierzenie się z różnego rodzaju problemami statystycznymi i aktuarialnymi. Poruszane problemy: Statystyka w badaniach sondażowych rynku, metody reprezentacyjne, Teoria i rachunek prawdopodobieństwa, statystyka opisowa, teoria estymacji, testowanie hipotez statystycznych, ekonometria, prognozowanie, metody data mining.
Copyright (C) 2006-2013 Statystycy.pl
Powered by phpBB modified by Przemo © 2003 phpBB Group
Strona wygenerowana w 0,12 sekundy. Zapytań do SQL: 11