În cartea mea recentă (vezi mai jos), la pagina 166 și mai devreme, am subliniat că, prin comparații pe perechi și, mai general, ori de câte ori se efectuează teste statistice simultane, este necesar să se furnizeze valori P care să țină seama de rata de eroare la nivel de familie, adică probabilitatea de a comite cel puțin o respingere incorectă în întreaga familie de valori P ajustate simultan. În acest sens, poate fi util să reamintim că, pentru un singur test nesemnificativ, rata de eroare din punct de vedere al comparației (E_c) este probabilitatea unei respingeri greșite pentru acel singur test (pe baza unei valori P neajustate), în timp ce probabilitatea a cel puțin unei respingeri greșite într-o familie de (k) comparațiile sunt mult mai mari.
Cu comparațiile perechi, un singur test se bazează de obicei pe raportul dintre o diferență și eroarea sa standard (un test t), care se presupune că urmează o distribuție t univariată atunci când ipoteza nulă este adevărată. Când sunt efectuate mai multe teste t simultane, se poate presupune că vectorul tuturor rapoartelor t urmează o distribuție t multivariată sub ipoteza că valoarea nulă este adevărată pentru toate testele simultane (Bretz et al., 2011). Prin urmare, valorile P ajustate pot fi obținute utilizând funcția de probabilitate a unei distribuții t multivariate în locul distribuției t simple univariate.
Ca exemplu, să reconsiderăm datele „amestec” utilizate în capitolul 9 al cărții principale. Trei amestecuri de erbicide și un martor netratat au fost testate pentru capacitatea lor de combatere a buruienilor împotriva unei buruieni importante din tomate, și anume Solanum nigrum. În codul de mai jos, încărcăm datele și adaptăm un model ANOVA unidirecțional, folosind greutatea plantelor de buruieni per ghiveci ca variabilă de răspuns și tratamentul cu erbicid ca factor explicativ. De dragul simplității, omitem verificările obișnuite ale ipotezelor de bază (vezi cartea principală). Tabelul ANOVA arată că efectul tratamentului este semnificativ și, prin urmare, trecem la compararea mijloacelor de tratament într-un mod perechi. Valorile P prezentate mai jos nu iau în considerare rata de eroare la nivel de familie, ci doar rata de eroare la nivel de comparație; aceste valori P pot fi reproduse folosind funcția de probabilitate a unei distribuții t Student univariate (pt() funcția în R).
library(statforbiology)
library(emmeans)
library(multcomp)
dataset <- getAgroData("mixture")
dataset$Treat <- factor(dataset$Treat)
model <- lm(Weight ~ Treat, data = dataset)
anova(model)
## Analysis of Variance Table
##
## Response: Weight
## Df Sum Sq Mean Sq F value Pr(>F)
## Treat 3 1089.53 363.18 23.663 2.509e-05 ***
## Residuals 12 184.18 15.35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
groupMeans <- emmeans(model, ~Treat)
tab <- contrast(groupMeans, method = "pairwise", adjust = "none")
tab
## contrast estimate SE df t.ratio p.value
## Metribuzin__348 - Mixture_378 4.05 2.77 12 1.461 0.1697
## Metribuzin__348 - Rimsulfuron_30 -7.68 2.77 12 -2.774 0.0168
## Metribuzin__348 - Unweeded -17.60 2.77 12 -6.352 <.0001
## Mixture_378 - Rimsulfuron_30 -11.73 2.77 12 -4.235 0.0012
## Mixture_378 - Unweeded -21.64 2.77 12 -7.813 <.0001
## Rimsulfuron_30 - Unweeded -9.91 2.77 12 -3.578 0.0038
#
# The P-value is obtained from the univariate t distribution (two-tails test)
abst <- abs(as.data.frame(tab)$t.ratio)
2 * pt(abst, 12, lower.tail = FALSE)
## (1) 1.696785e-01 1.683167e-02 3.651239e-05 1.157189e-03 4.782986e-06
## (6) 3.794451e-03
Pentru a obține rate de eroare la nivel de familie, ar trebui să trecem de la distribuția t univariată la cea multivariată. De exemplu, să luăm în considerare primul raport t din caseta de cod anterioară (t = 1,461). Ar trebui să ne întrebăm: „care este probabilitatea de a obține un raport t la fel de extrem sau mai extrem decât 1,461 dintr-o distribuție t multivariată cu șase dimensiuni (adică numărul de teste simultane)?”. În acest calcul, mai trebuie să avem în vedere că cele 6 teste sunt corelate, cel puțin într-o oarecare măsură, pentru că au câteva elemente comune, de exemplu, același termen de eroare la numitor. În cel mai simplu caz (homoscedasticitate și date echilibrate), această corelație este egală cu 0,5 pentru toate comparațiile pe perechi.
În vremuri mai vechi, când puterea de calcul era limitată, calcularea probabilităților din distribuția t multivariată era o sarcină descurajantă. Cu toate acestea, pentru unele cazuri specifice (de exemplu, modele liniare cu date homoscedastice și echilibrate), valorile P ajustate pot fi obținute prin exploatarea distribuției intervalului studentizat (așa-numita metodă „tukey”), care este opțiunea implicită în contrast() funcția de emmeans pachet, așa cum se arată în caseta de cod următoare.
tab <- contrast(groupMeans, method = "pairwise") # tab <- contrast(groupMeans, method = "pairwise", adjust = "tukey") # same as above tab ## contrast estimate SE df t.ratio p.value ## Metribuzin__348 - Mixture_378 4.05 2.77 12 1.461 0.4885 ## Metribuzin__348 - Rimsulfuron_30 -7.68 2.77 12 -2.774 0.0698 ## Metribuzin__348 - Unweeded -17.60 2.77 12 -6.352 0.0002 ## Mixture_378 - Rimsulfuron_30 -11.73 2.77 12 -4.235 0.0055 ## Mixture_378 - Unweeded -21.64 2.77 12 -7.813 <.0001 ## Rimsulfuron_30 - Unweeded -9.91 2.77 12 -3.578 0.0173 ## ## P value adjustment: tukey method for comparing a family of 4 estimates # The P-value is obtained from the Studentised Range Distribution (two-tails test) abst <- abs(as.data.frame(tab)$t.ratio) ptukey(sqrt(2) * abst, 4, 12, lower.tail = FALSE) ## (1) 4.884620e-01 6.981178e-02 1.853807e-04 5.501451e-03 2.473776e-05 ## (6) 1.725725e-02
Această metodă simplă oferă rate exacte de eroare la nivel de familie cu date echilibrate – care reprezintă marea majoritate a experimentelor de teren concepute în agricultură – și funcționează rezonabil de bine în prezența unor grade mici de dezechilibru. În cadrul procedurilor tradiționale de testare cu comparații multiple, abordarea descrisă mai sus conduce la aceleași rezultate ca și HSD-ul Tukey pentru date echilibrate și testul Tukey-Kramer pentru date neechilibrate.
Mai recent, a devenit posibil să se calculeze direct probabilitățile din distribuția t multivariată, ceea ce este deosebit de convenabil deoarece oferă o abordare mai generală pentru obținerea ratelor de eroare la nivel de familie. Această distribuție este implementată în pachetul „mvtnorm” prin intermediul pmvt() funcţie. Pentru a efectua calculul, trebuie să precizăm, pentru fiecare dimensiune, intervalul peste care urmează să fie calculată probabilitatea (în acest caz, pentru primul raport t, intervalul este (pm 1,461081)), numărul de grade de libertate (12) și matricea de corelație a combinațiilor liniare, care pot fi preluate direct din obiectul „emmGrid”. Codul de mai jos ilustrează aceste calcule. Mărimea „plev” reprezintă probabilitatea de eșantionare în cadrul intervalului (adică nici una dintre cele șase ipoteze nule nu este respinsă greșit), în timp ce rata de eroare la nivel de familie corespunde probabilității de eșantionare în afara intervalului (adică cel puțin o ipoteză nulă este respinsă greșit), care se obține prin scădere.
library(mvtnorm)
t1 <- abs(as.data.frame(tab)$t.ratio)(1)
ncontr <- 6
corMat <- cov2cor(vcov(tab))
plev <- pmvt(lower = rep(-t1, ncontr), upper=rep(t1, ncontr), df = 12,
corr = corMat)(1)
1 - plev
## (1) 0.4883843
În R, o astfel de abordare poate fi obținută utilizând adjust = "mvt" argument.
tab <- contrast(groupMeans, method = "pairwise", adjust = "mvt") tab ## contrast estimate SE df t.ratio p.value ## Metribuzin__348 - Mixture_378 4.05 2.77 12 1.461 0.4885 ## Metribuzin__348 - Rimsulfuron_30 -7.68 2.77 12 -2.774 0.0698 ## Metribuzin__348 - Unweeded -17.60 2.77 12 -6.352 0.0002 ## Mixture_378 - Rimsulfuron_30 -11.73 2.77 12 -4.235 0.0054 ## Mixture_378 - Unweeded -21.64 2.77 12 -7.813 <.0001 ## Rimsulfuron_30 - Unweeded -9.91 2.77 12 -3.578 0.0172 ## ## P value adjustment: mvt method for 6 tests
Funcția de mai sus folosește metode de integrare numerică și se bazează pe simulare; în consecință, rezultatele nu sunt pe deplin reproductibile. Cu toate acestea, este ușor de observat că aceste rezultate sunt echivalente asimptotic cu cele obținute cu metoda de ajustare Tukey prezentată mai sus. Datorită acestei complexități intrinseci, utilizarea adjust = "mvt" argumentul nu este recomandat pentru comparații în perechi în experimente echilibrate, în timp ce se poate dovedi util în alte situații, de exemplu în prezența unor date puternic dezechilibrate.
Vă mulțumim pentru citit și nu uitați să vedeți noua mea carte de mai jos!
Andrea Onofri
Departamentul de Științe Agricole, Alimentației și Mediului
Universitatea din Perugia (Italia)
Trimite comentarii la: (e-mail protejat)
