Într-o postare recentă am arătat că putem construi combinații liniare ale parametrilor modelului (vezi aici ). De exemplu, dacă avem două estimări ale parametrilor, să spunem Q și W, cu erori standard, respectiv egale cu (sigma_Q) şi (sigma_W)putem construi o combinație liniară după cum urmează:
( Z = aQ + bW + c)
unde (o), (b) şi (c) sunt trei coeficienți. Eroarea standard pentru această combinație poate fi obținută ca:
( sigma_Z = sqrt{ a^2 sigma^2_Q + b^2 sigma^2_W + 2ab sigma_{QW} })
În biologie, transformările neliniare sunt mult mai frecvente decât transformările liniare. Transformările neliniare sunt, de exemplu, (Z = exp(Q + W))sau, (Z = 1/(Q + W)). Care este eroarea standard pentru aceste transformări neliniare? Aceasta nu este o problemă complexă, dar soluția poate fi dincolo de biologii cu un nivel mediu de competență statistică. Este numită „metoda delta” și oferă așa-numitele „erori standard delta”. M-am gândit că ar putea fi util să vorbesc despre asta, folosind un limbaj foarte simplu și câteva exemple.
Un erbicid s-a dovedit a urma o cinetică de degradare de prim ordin în sol, cu o rată constantă de degradare (k = -0,035) și eroare standard egală cu (0,00195). Care este timpul de înjumătățire ((T_{1/2})) a acestui erbicid și eroarea sa standard?
Fiecare chimist de pesticide știe că timpul de înjumătățire ((T_{1/2})) este derivată din rata de degradare, conform următoarei ecuații:
( T_{1/2} = frac{log(0,5)}{k})
Prin urmare, timpul de înjumătățire al erbicidului nostru este:
Thalf <- log(0.5)/-0.035 Thalf ## (1) 19.80421
Dar… care este eroarea standard a acestui timp de înjumătățire? Există o oarecare incertitudine în jurul estimării (k) şi este clar că o astfel de incertitudine ar trebui să se propage la estimarea a (T_{1/2}); din păcate, transformarea este neliniară și nu putem folosi expresia dată mai sus pentru transformări liniare.
Ideea de bază
Ideea de bază din spatele metodei deltei este că majoritatea funcțiilor neliniare simple, pe care le folosim în biologie, pot fi aproximate local prin linia tangentă printr-un punct de interes. De exemplu, funcția noastră neliniară de înjumătățire este (Y = frac{log(0,5)}{X}) și, evident, ne interesează punctul în care (X = k = -0,035). În graficul de mai jos, am reprezentat funcția noastră neliniară (cu negru) și linia ei tangentă (cu roșu) prin punctul de mai sus: putem vedea că aproximarea este destul de bună în imediata apropiere a (X = -0,035).
Care este ecuația dreptei tangente? În general, dacă funcția neliniară este (G(X))poate vă amintiți din liceu că pârtia (m) a dreptei tangente este egală cu prima derivată a (G(X))adică (G'(X)). Vă puteți aminti, de asemenea, că ecuația unei drepte cu panta (m) prin punct (P(X_1, Y_1)) este (Y – Y_1 = m(X – X_1)). Ca (Y_1 = G(X_1))linia tangentă are ecuația:
(Y = G(X_1) + G'(X_1)(X – X_1))
Avem nevoie de derivat!
Pentru a scrie ecuația liniei roșii din figura de mai sus, trebuie să luăm în considerare asta (X_1 = -0,035) și trebuie să putem calcula prima derivată a funcției noastre de înjumătățire neliniară. Nu pot deriva expresia primei derivate pentru toate funcțiile neliniare și a fost o ușurare pentru mine să descopăr că R poate face față acestei sarcini în moduri simple, de exemplu folosind funcția D(). Pentru cazul nostru, este:
D(expression(log(0.5)/X), "X") ## -(log(0.5)/X^2)
Prin urmare, putem folosi această funcție R pentru a calcula panta (m) a dreptei tangente din figura de mai sus:
X <- -0.035 m <- eval( D(expression(log(0.5)/X), "X") ) m ## (1) 565.8344
Știm deja asta (G(-0,035) = 19,80421). Prin urmare, putem scrie ecuația dreptei tangente (linia roșie în graficul de mai sus):
(Y = 19,80421 + 565,8344 , (X + 0,035))
adica:
(Y = 19,80421 + 565,8344 , X + 565,8344 cdot 0,035 = 39,60841 + 565,8344 , X)
Înlocuirea unei curbe cu o linie
Acum avem două funcții:
- funcția de înjumătățire neliniară inițială (Y = log(0,5)/X)$
- o nouă funcție liniară ((Y = 39,60841 + 565,8344 , X)), adică o aproximare foarte apropiată de cea precedentă, cel puțin aproape de punct (X = -0,035)care ne interesează.
Prin urmare, îl putem aproxima pe primul cu cel din urmă! Dacă folosim funcția liniară, vedem că timpul de înjumătățire este:
39.60841 + 565.8344 * -0.035 ## (1) 19.80421
care este ceea ce ne așteptam. Avantajul este că acum putem folosi valoarea scăzută de propagare a erorilor pentru a estima eroarea standard (vezi prima și a doua ecuație din această postare):
( sigma_{ left( 39,60841 + 565,8344 , X right)} = sqrt{ 562,8344^2 , sigma^2_X})
Începem:
sqrt( m^2 * (0.00195 ^ 2) ) ## (1) 1.103377
In general…
Dacă avem o transformare neliniară (G(X))eroarea standard pentru această transformare este aproximată prin cunoașterea primei derivate (G'(X)) iar eroarea standard a (X):
(sigma_{G(X)} simeq sqrt{ (G'(X))^2 , sigma^2_X })
O lucrare raportează că numărul mediu de microorganisme dintr-un substrat, la scară logaritmică, a fost (X_1 = 5) cu eroare standard (sigma = 0,84). Este ușor de dedus că numărul real de microorganisme a fost (exp{5} = 148,4132); care este eroarea standard a mediei transformate înapoi?
Prima derivată a funcției noastre neliniare este:
D(expression(exp(X)), "X") ## exp(X)
și astfel panta dreptei tangente este:
X <- 5 m <- eval( D(expression(exp(X)), "X") ) m ## (1) 148.4132
Conform funcției de mai sus, eroarea standard pentru media transformată înapoi este:
sigma <- 0.84 sqrt( m^2 * sigma^2 ) ## (1) 124.6671
Concentrația de seleniu în drupele de măsline s-a dovedit a fi (3,1 , mu g ,, g^{-1}) cu eroare standard egală cu 0,8. Care este aportul de seleniu atunci când mănânci o drupă? Vă rugăm să luați în considerare că o drupă cântărește, în medie, 3,4 g (SE = 0,31) și că concentrația de seleniu și greutatea drupei prezintă o covarianță de 0,55.
Cantitatea de seleniu este ușor de calculată astfel:
X <- 3.1; W = 3.4 X * W ## (1) 10.54
Erorile standard Delta pot fi obținute luând în considerare derivatele parțiale pentru fiecare dintre cele două variabile:
mX <- eval( D(expression(X*W), "X") ) mW <- eval( D(expression(X*W), "W") )
și combinându-le după cum urmează:
sigmaX <- 0.8; sigmaW <- 0.31; sigmaXW <- 0.55 sqrt( (mX^2) * sigmaX^2 + (mW^2) * sigmaW^2 + 2 * X * W * sigmaXW ) ## (1) 4.462726
Pentru aceia dintre voi care ar dori să se implice în notația matriceală: putem ajunge la același rezultat prin multiplicarea matricei (vezi mai jos). Acest lucru ar putea fi mai ușor atunci când avem mai mult de două variabile de combinat.
der <- matrix(c(mX, mW), 1, 2) sigma <- matrix(c(sigmaX^2, sigmaXW, sigmaXW, sigmaW^2), 2, 2, byrow = T) sqrt( der %*% sigma %*% t(der) ) ## (,1) ## (1,) 4.462726
În R există o funcție de comandă rapidă pentru a calcula erorile standard delta, care este disponibilă în pachetul „mașină”. Pentru a-l folosi, trebuie să avem:
- un vector numit pentru variabilele pe care trebuie să le combinăm
- o expresie pentru transformare
- o matrice varianță-covarianță
Pentru primul exemplu avem:
obj <- c("k" = -0.035)
sigma <- matrix(c(0.00195^2), 1, 1)
library(car)
## Loading required package: carData
deltaMethod(object = obj, g="log(0.5)/k", vcov = sigma)
## Estimate SE 2.5 % 97.5 %
## log(0.5)/k 19.8042 1.1034 17.6416 21.967
Pentru al doilea exemplu:
obj <- c("X1" = 5)
sigma <- matrix(c(0.84^2), 1, 1)
deltaMethod(object = obj, g="exp(X1)", vcov = sigma)
## Estimate SE 2.5 % 97.5 %
## exp(X1) 148.41 124.67 -95.93 392.76
Pentru al treilea exemplu:
obj <- c("X" = 3.1, "W" = 3.4)
sigma <- matrix(c(0.8^2, 0.55, 0.55, 0.31^2), 2, 2, byrow = T)
deltaMethod(object = obj, g="X * W", vcov = sigma)
## Estimate SE 2.5 % 97.5 %
## X * W 10.5400 4.4627 1.7932 19.287
Funcția „deltaMethod()” este foarte utilă pentru a fi folosită în legătură cu obiectele model, deoarece nu trebuie să furnizăm nimic, decât funcția de transformare. Dar acesta este ceva care necesită o altă postare!
Cu toate acestea, două note finale referitoare la metoda delta trebuie subliniate aici:
- eroarea standard delta este întotdeauna aproximativă;
- dacă variabilele originale sunt gaussiene, variabila transformată, de obicei, nu este gaussiană.
Multumesc pentru lectura!
prof. Andrea Onofri
Departamentul de Științe Agricole, Alimentației și Mediului
Universitatea din Perugia (Italia)
Trimite comentarii la: (e-mail protejat)
