Dacă observațiile ar fi independente (adică fără blocuri și fără măsuri repetate), acest model ar putea fi potrivit utilizând regresia neliniară convențională. Preferința mea merge la drm() funcția în pachetul „drc” (Ritz și colab., 2015).
Codarea este raportată mai jos: „Randamentul” este o funcție a ((sim)) DAS, prin intermediul unei funcții logistice cu trei parametri („fct = L.3()”). Trebuie ajustate curbe diferite pentru diferite combinații de genotip și niveluri de azot („curbeid = N:GEN”), deși aceste curbe ar trebui să se bazeze parțial pe valorile parametrilor comuni („pmodele = …).
Argumentul „pmodels” necesită câteva comentarii suplimentare. Trebuie să fie un vector cu atâtea elemente câte parametri există în model (trei, în acest caz: (b), (d) şi (e)). Fiecare element reprezintă o funcție liniară a variabilelor și se referă la parametrii în ordine alfabetică, adică primul element se referă la (b)al doilea se referă la (d) iar al treilea să (e). Parametrul (b) nu este dependent de nicio variabilă (‘~ 1’) și astfel o valoare constantă este ajustată peste curbe; (d) şi (e) depind de o combinație complet factorială de genotip și nivel de azot (~ N*GEN = ~N + GEN + N:GEN). În cele din urmă, am folosit argumentul „bcVal = 0,5” pentru a specifica că intenționăm să folosim o tehnică Transform-Both-Sides, în care se efectuează o transformare logaritmică pentru ambele părți ale ecuației. Acest lucru este necesar pentru a lua în considerare heteroscedasticitatea, dar nu afectează scara estimărilor parametrilor.
Acest model poate fi util pentru alte circumstanțe (fără blocuri și fără măsuri repetate), dar este greșit în exemplul nostru. Într-adevăr, observațiile sunt grupate în blocuri și parcele; neglijând acest lucru, încălcăm ipoteza independenței reziduurilor modelului. Graficele rapide ale reziduurilor față de valorile ajustate și diagrama QQ arată că nu există probleme cu heteroscedasticitatea și normalitatea, deși câteva valori aberante ar putea merita o inspecție atentă (pe care le vom neglija acum). Rețineți că aceste două grafice sunt obținute folosind plotRes() funcția în „statforbiologie”.
Având în vedere cele de mai sus, trebuie să folosim un alt model, aici, deși voi arăta că această potrivire naivă se poate dovedi utilă.
Potrivire de model mixt neliniar
Pentru a ține seama de gruparea observațiilor, trecem la un model neliniar cu efect mixt (NLME). O alegere bună este nlme() funcția în pachetul „nlme” (Pinheiro și Bates, 2000), deși sintaxa poate fi greoaie, uneori. Voi încerca să ajut, enumerând și comentând cele mai importante argumente pentru această funcție. Trebuie să precizăm următoarele:
- O funcție deterministă. În acest caz, folosim
NLS.L3()funcția în pachetul „statforbiology”, care oferă un model de creștere logistică cu aceeași parametrizare ca șiL.3()funcția din pachetul „drc”. - Funcții liniare pentru parametrii modelului. Argumentul „fix” din funcția „nlme” este foarte asemănător cu argumentul „pmodels” din funcția „drm” de mai sus, în sensul că necesită o listă, în care fiecare element este o funcție liniară a variabilelor. Singura diferență este că numele parametrului trebuie specificat în partea stângă a funcției.
- Efecte aleatorii pentru parametrii modelului. Acestea sunt specificate folosind argumentul „aleatoriu”. În acest caz, parametrii (d) şi (e) se așteaptă să prezinte variabilitate aleatorie de la un bloc la altul și de la un teren la altul, în cadrul unui bloc. De dragul simplității, ca parametru (b) nu este afectat de genotip și nivelul de azot, ne așteptăm, de asemenea, să nu prezinte nicio variabilitate aleatorie între blocuri și parcele.
- Valori de pornire pentru parametrii modelului. Rutinele de autopornire nu sunt folosite de
nlme()și astfel trebuie să specificăm un vector numit, care să păstreze valorile inițiale ale parametrilor modelului. În acest caz, am decis să folosesc rezultatul din regresia neliniară „naivă” de mai sus, care, prin urmare, se dovedește utilă.
Transformarea ambelor părți ale ecuației se face în mod explicit.
modnlme1 <- nlme(sqrt(Yield) ~ sqrt(NLS.L3(DAS, b, d, e)),
data = dataset,
random = d + e ~ 1|Block/Plot,
fixed = list(b ~ 1, d ~ N*GEN, e ~ N*GEN),
start = coef(modNaive1),
control = list(msMaxIter = 400))
summary(modnlme1)$tTable
## Value Std.Error DF t-value p-value
## b 0.05652837 0.002157629 228 26.1993018 1.044743e-70
## d.(Intercept) 33.91575348 1.222612866 228 27.7403865 5.073982e-75
## d.NLow -3.18382830 1.592502935 228 -1.9992606 4.676721e-02
## d.GENB 18.90014651 1.864712716 228 10.1356881 3.602004e-20
## d.GENC -1.15934033 1.686405282 228 -0.6874625 4.924901e-01
## d.NLow:GENB -5.99216682 2.455759119 228 -2.4400467 1.544863e-02
## d.NLow:GENC -5.82864845 2.217534809 228 -2.6284361 9.160826e-03
## e.(Intercept) 55.20071252 2.320784607 228 23.7853665 1.087153e-63
## e.NLow -9.06217426 3.127165879 228 -2.8978873 4.123120e-03
## e.GENB -4.47038042 2.761533544 228 -1.6188036 1.068720e-01
## e.GENC 4.00746125 3.084383627 228 1.2992746 1.951620e-01
## e.NLow:GENB -4.71367457 4.055953633 228 -1.1621618 2.463848e-01
## e.NLow:GENC 2.23951067 4.609547644 228 0.4858417 6.275460e-01
Din diagramele de mai sus, vedem că potrivirea generală este bună. Efectele fixe și componentele de varianță pentru efecte aleatoare sunt obținute după cum urmează:
summary(modnlme1)$tTable ## Value Std.Error DF t-value p-value ## b 0.05652837 0.002157629 228 26.1993018 1.044743e-70 ## d.(Intercept) 33.91575348 1.222612866 228 27.7403865 5.073982e-75 ## d.NLow -3.18382830 1.592502935 228 -1.9992606 4.676721e-02 ## d.GENB 18.90014651 1.864712716 228 10.1356881 3.602004e-20 ## d.GENC -1.15934033 1.686405282 228 -0.6874625 4.924901e-01 ## d.NLow:GENB -5.99216682 2.455759119 228 -2.4400467 1.544863e-02 ## d.NLow:GENC -5.82864845 2.217534809 228 -2.6284361 9.160826e-03 ## e.(Intercept) 55.20071252 2.320784607 228 23.7853665 1.087153e-63 ## e.NLow -9.06217426 3.127165879 228 -2.8978873 4.123120e-03 ## e.GENB -4.47038042 2.761533544 228 -1.6188036 1.068720e-01 ## e.GENC 4.00746125 3.084383627 228 1.2992746 1.951620e-01 ## e.NLow:GENB -4.71367457 4.055953633 228 -1.1621618 2.463848e-01 ## e.NLow:GENC 2.23951067 4.609547644 228 0.4858417 6.275460e-01 # VarCorr(modnlme1) ## Variance StdDev Corr ## Block = pdLogChol(list(d ~ 1,e ~ 1)) ## d.(Intercept) 3.933210e-08 1.983232e-04 d.(In) ## e.(Intercept) 1.747130e-08 1.321791e-04 -0.001 ## Plot = pdLogChol(list(d ~ 1,e ~ 1)) ## d.(Intercept) 3.197539e-09 5.654679e-05 d.(In) ## e.(Intercept) 9.573313e-10 3.094077e-05 0 ## Residual 1.750623e-01 4.184044e-01
Acum, să ne întoarcem la scopul nostru inițial: testarea semnificației interacțiunii „genotip x azot”. Într-adevăr, avem două teste disponibile: pornit pentru parametru (d) și unul pentru parametru (e). În primul rând, codificăm două modele „reduse”, în care efectele genotipului și azotului sunt pur dependente. Pentru a face acest lucru, schimbăm specificația efectelor fixe de la „~ N*GEN” la „~ N + GEN”. De asemenea, în acest caz, folosim o potrivire de regresie neliniară „naivă” pentru a obține valori inițiale pentru parametrii modelului, care vor fi utilizate în următoarea potrivire a modelului NLME.
modNaive2 <- drm(Yield ~ DAS, fct = L.3(),
data = dataset,
curveid = N:GEN,
pmodels = c( ~ 1, ~ N + GEN, ~ N * GEN),
bcVal = 0.5)
modnlme2 <- nlme(sqrt(Yield) ~ sqrt(NLS.L3(DAS, b, d, e)),
data = dataset,
random = d + e ~ 1|Block/Plot,
fixed = list(b ~ 1, d ~ N + GEN, e ~ N*GEN),
start = coef(modNaive2),
control = list(msMaxIter = 200))
modNaive3 <- drm(Yield ~ DAS, fct = L.3(), data = dataset,
curveid = N:GEN,
pmodels = c( ~ 1, ~ N*GEN, ~ N + GEN), bcVal = 0.5)
modnlme3 <- nlme(sqrt(Yield) ~ sqrt(NLS.L3(DAS, b, d, e)),
data = dataset,
random = d + e ~ 1|Block/Plot,
fixed = list(b ~ 1, d ~ N*GEN, e ~ N + GEN),
start = coef(modNaive3), control = list(msMaxIter = 200))
Să luăm în considerare primul model redus „modnlme2”. În acest model, interacțiunea „genotip x azot” a fost eliminată pentru (d) parametru. Putem compara acest model redus cu modelul complet „modnlme1”, folosind un test al raportului de probabilitate:
anova(modnlme1, modnlme2) ## Model df AIC BIC logLik Test L.Ratio p-value ## modnlme1 1 20 329.1496 400.6686 -144.5748 ## modnlme2 2 18 334.2187 398.5857 -149.1093 1 vs 2 9.069075 0.0107
Acest test este semnificativ, dar valoarea AIC este foarte apropiată pentru cele două modele. Având în vedere că un LRT în modele mixte este de obicei destul de liberal, ar trebui să se poată concluziona că interacțiunea „genotip x azot” nu este semnificativă și, prin urmare, clasarea genotipurilor în ceea ce privește potențialul de randament, măsurată prin (d) parametrul ar trebui să fie independent de nivelul de azot.
Să luăm acum în considerare al doilea model redus „modnlme3”. În acest al doilea model, interacțiunea „genotip x azot” a fost eliminată pentru parametrul „e”. Putem compara și acest model redus cu modelul complet „modnlme1”, utilizând un test al raportului de probabilitate:
anova(modnlme1, modnlme3) ## Model df AIC BIC logLik Test L.Ratio p-value ## modnlme1 1 20 329.1496 400.6686 -144.5748 ## modnlme3 2 18 328.2446 392.6117 -146.1223 1 vs 2 3.095068 0.2128
În acest al doilea test, lipsa de semnificație pentru interacțiunea „genotip x azot” pare să fie mai puțin discutabilă decât în primul.
Aș dori să închei prin a vă atrage atenția asupra funcției „medrm” din pachetul „medrc”, care poate fi folosită și pentru a se potrivi cu acest tip de modele neliniare cu efecte mixte.
Codare fericită cu R! … și … nu uitați să verificați noua mea carte!
prof. Andrea Onofri
Departamentul de Științe Agricole, Alimentației și Mediului
Universitatea din Perugia (Italia)
Trimite comentarii la: (e-mail protejat)
