În activitatea noastră de cercetare, de obicei potrivim modele la datele experimentale. Scopul nostru este de a estima unii parametri biologic relevanți, împreună cu erorile standard ale acestora. Foarte des, acești parametri sunt interesanți în sine, deoarece reprezintă medii, diferențe, rate sau alți descriptori importanți. În alte cazuri, folosim acele estimări pentru a obține indici suplimentari, prin intermediul unor calcule adecvate. De exemplu, gândiți-vă că avem două estimări ale parametrilor, să spunem Q și W, cu erori standard, respectiv egale cu (sigma_Q) şi (sigma_W): ar putea fi relevant să se calculeze suma:
(Z = aQ + bW + c)
unde (o), (b) şi (c) sunt trei coeficienți. Operația de mai sus se numește „combinație liniară”; este un fel de sumă ponderată a parametrilor modelului. Întrebarea este: care este eroarea standard a lui Z? Aș dori să arăt acest lucru printr-un exemplu biologic simplu.
Am măsurat rata de germinare a semințelor de Brassica rapa la șase nivele de potențial de apă în substrat (-1, -0,9, -0,8, -0,7, -0,6 și -0,5 MPa). Am dori să prezicăm rata de germinare pentru un nivel de potențial de apă de -0,75 MPa.
Referințele din literatură sugerează că relația dintre rata de germinare și potențialul de apă din substrat este liniară. Prin urmare, ca prim pas, adaptăm un model de regresie liniară la datele noastre observate. Dacă suntem în R, codul pentru a îndeplini această sarcină este afișat mai jos.
GR <- c(0.11, 0.20, 0.34, 0.46, 0.56, 0.68) Psi <- c(-1, -0.9, -0.8, -0.7, -0.6, -0.5) lMod <- lm(GR ~ Psi) summary(lMod) ## ## Call: ## lm(formula = GR ~ Psi) ## ## Residuals: ## 1 2 3 4 5 6 ## 0.0076190 -0.0180952 0.0061905 0.0104762 -0.0052381 -0.0009524 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.25952 0.02179 57.79 5.37e-07 *** ## Psi 1.15714 0.02833 40.84 2.15e-06 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.01185 on 4 degrees of freedom ## Multiple R-squared: 0.9976, Adjusted R-squared: 0.997 ## F-statistic: 1668 on 1 and 4 DF, p-value: 2.148e-06
Este clar că putem folosi modelul adaptat pentru a calcula valoarea GR la -0,75 MPa, ca (GR = 1,26 – 1,16 time 0,75 = 0,39). Aceasta este într-adevăr o combinație liniară a parametrilor modelului, așa cum am arătat mai sus. Q și W sunt parametrii modelului estimați, în timp ce (a = 1), (b = -0,75) şi (c = 0).
Desigur, răspunsul derivat este, de asemenea, o estimare și trebuie să dăm o măsură a incertitudinii. Pentru ambii parametri de model avem erori standard; întrebarea este: cum se propagă incertitudinea în estimările parametrilor la combinația lor liniară? Cu cuvinte mai simple: este ușor să construim o sumă ponderată a parametrilor modelului, dar cum facem o sumă ponderată a erorilor standard ale acestora?
Sokal și Rohlf (1981) la pag. 818 din cartea lor, arată că, în general, este:
( textrm{var}(a , Q + b , W) = a^2 sigma^2_Q + b^2 sigma^2_W + 2ab sigma_{QW} )
unde (sigma_{QW}) este covarianța lui Q și W. Atașez mai jos demonstrația; este destul de simplu de înțeles și merită efortul. Cu toate acestea, mai mulți studenți la biologie sunt destul de reticenți atunci când trebuie să se aprofundeze în matematică. Prin urmare, aș dori, de asemenea, să ofer o „dovadă” empirică, folosind un cod R simplu.
Să luăm două mostre (Q și W) și să le combinăm folosind doi coeficienți (o) şi (b). Să calculăm varianța pentru combinație:
Q <- c(12, 14, 11, 9) W <- c(2, 4, 7, 8) a <- 2 b <- 3 c <- 4 var(a * Q + b * W + c) ## (1) 35.58333 r a^2 * var(Q) + b^2 * var(W) + 2 * a * b * cov(Q, W) ## (1) 35.58333
Vedem ca rezultatele se potrivesc exact! În exemplul nostru, varianța și covarianța estimărilor sunt preluate folosind funcția „vcov()”:
vcov(lMod) ## (Intercept) Psi ## (Intercept) 0.0004749433 0.0006020408 ## Psi 0.0006020408 0.0008027211 r sigma2Q <- vcov(lMod)(1,1) sigma2W <- vcov(lMod)(2,2) sigmaQW <- vcov(lMod)(1,2)
Eroarea standard pentru predicție se obține pur și simplu ca:
sqrt( sigma2Q + 0.75^2 * sigma2W - 2 * 0.75 * sigmaQW ) ## (1) 0.004838667
Acum că am înțeles conceptul, putem folosi R pentru a face calculele. Pentru acest exemplu, metoda „predict()” reprezintă cea mai logică abordare. Trebuie doar să trecem obiectul model și valoarea X pentru care trebuie să facem o predicție. Această din urmă valoare trebuie organizată ca un cadru de date, cu numele coloanelor care să se potrivească cu numele variabilelor X din setul de date original:
predict(lMod, newdata = data.frame(Psi = -0.75),
se.fit = T)
## $fit
## 1
## 0.3916667
##
## $se.fit
## (1) 0.004838667
##
## $df
## (1) 4
##
## $residual.scale
## (1) 0.01185227
În afară de metoda predicției, există o altă funcție de interes mai general, care poate fi folosită pentru a construi combinații liniare ale parametrilor modelului. Este funcția „glht()” din pachetul „multcomp”. Pentru a folosi această funcție avem nevoie de un obiect model și trebuie să organizăm coeficienții pentru transformare într-o matrice, cu atâtea rânduri câte combinații sunt de calculat. Când scriem coeficienții, ignorăm C: am văzut că fiecare valoare constantă nu afectează varianța transformării.
De exemplu, imaginați-vă că vrem să prezicăm GR pentru două niveluri de potențial de apă, adică -0,75 (ca mai sus) și -0,55 MPa. Coeficienții (A, B) pentru combinații sunt:
Z1 <- c(1, -0.75) Z2 <- c(1, -0.55)
Adunăm cei doi vectori într-o singură matrice cu două rânduri:
M <- matrix(c(Z1, Z2), 2, 2, byrow = T)
Și acum trecem acea matrice funcției „glht()” ca argument:
library(multcomp) lcombs <- glht(lMod, linfct = M) summary(lcombs) ## ## Simultaneous Tests for General Linear Hypotheses ## ## Fit: lm(formula = GR ~ Psi) ## ## Linear Hypotheses: ## Estimate Std. Error t value Pr(>|t|) ## 1 == 0 0.391667 0.004839 80.94 2.30e-07 *** ## 2 == 0 0.623095 0.007451 83.62 2.02e-07 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## (Adjusted p values reported -- single-step method)
Metoda de mai sus poate fi extinsă cu ușurință la alte tipuri de modele liniare și combinații liniare de parametri de model. Argumentul „ajustare” este util ori de câte ori dorim să obținem intervale de încredere la nivel de familie, care pot explica problema multiplicității. Dar asta e alta poveste…
Lucru fericit cu aceste funcții!
prof. Andrea Onofri
Departamentul de Științe Agricole, Alimentației și Mediului
Universitatea din Perugia (Italia)
Trimite comentarii la: (e-mail protejat)
Știm că varianța pentru o variabilă aleatoare este definită ca:
({begin{array}{l} var(Z) = frac{1}{n-1}sum left(Z – hat{Z}right)^2 = \ = frac{ 1}{n-1}sum left(Z – frac{1}{n} sum{Z}right)^2 end{array}})
Dacă (Z = aQ + bW + c)unde (o), (b) şi (c) sunt coeficienții și Q și W sunt două variabile aleatoare, avem:
( {begin{array}{l} var(Z) = frac{1}{n-1}sum left(aQ + bW + c – frac{1}{n} sum{ left (aQ + bW + cright)}right)^2 = \ = frac{1}{n-1}sum left(aQ + bW + c – frac{1}{n} sum {aQ} – frac{1}{n} sum{ bW} – frac{1}{n} sum{ c}right)^2 = end{array}} )
({begin{array}{l} = frac{1}{n-1}sum left(aQ + bW + c – a hat{Q} – b hat{W} – c right )^2 =\ = frac{1}{n-1}sum left(left( aQ – a hat{Q}right) + left( aW – ahat{W}right) ) right)^2 = end{matrice}})
({begin{array}{l} =frac{1}{n-1}sum left(a left( Q – hat{Q}right) + b left( W – hat) {W}right) right)^2 =\ = frac{1}{n-1}sum left(a^2 left( Q – hat{Q}right)^2 + b ^2 left( W – hat{W}right)^2 + 2ab left( Q – hat{Q}right) left( W – hat{W}right)right) = end{array}} )
( = a^2 frac{1}{n-1} sum{left( Q – hat{Q}right)^2} + b^2 frac{1}{n-1} suma left( W – hat{W}right)^2 + 2ab frac{1}{n-1}sum{left(left( Q – hat{Q}right) left( W – hat{W}dreapta)dreapta)})
Prin urmare:
( textrm{var}(Z) = a^2 sigma^2_Q + b^2 sigma^2_W + 2ab sigma_{Q,W})
