A lineáris modellektől a nemlineáris kevert modellekig R-ben

A lineáris modellektől a nemlineáris kevert modellekig R-ben Harnos Andrea Szent István Egyetem, Állatorvostudományi Kar Biomatematika Tanszék

Tartalom Bevezetés Modellezés Az általános lineáris modell Általánosított lineáris modellek Additív modellek Kevert modellek Nemlineáris kevert modell

Sweave

I

A kurzus diái Fritz Leisch’s Sweave rendszerével készültek. Ebben a LATEX és R kódok egyetlen fájlban szerkeszthetők.

I

Egy Sweave fájlból (.Rnw kiterjesztésű általában) olyan LATEX forrás fájl készíthető, amely tartalmazza az R inputokat, outputokat és ábrákat.

I

Egy Sweave fájlból az R kódok automatikusan kinyerhetők.

A kurzus anyagához felhasznált könyvek és egyéb anyagok

I

Brian S. Everitt, Torsten Hothorn: A Handbook of Statistical Analysis Using R (Chapman and Hall/CRC, 2006)

I

José C. Pinheiro, Douglas M. Bates: Mixed Effects Models in S and S-PLUS (Springer, 2000)

I

Julian J. Faraway: Extending the Linear Model with R. Generalized Linear, Mixed Effects and Nonparametric Regression Models (Chapman and Hall/CRC, 2006)

I

Douglas M. Bates: Mixed-effects models in R. useR!2006, Vienna, Austria, June 14, 2006

Felhasznált adatok

I

pupa.txt

I

mass.csv

I

land-use.csv

Zerynthia polyxena

pupa.txt, mass.csv I

Az adatok egy olyan kísérletből származnak, amelyben lepkék (Zerynthia polyxena) imágóinak méret változatosságát vizsgálták.

I

A lárvákat kísérletileg manipulált hőmérsékletű környezetben tartották.

I

A KÍSÉRLET 1. faktor: TEMPR Fejlődő hernyók környezetének hőmérséklete I I I

hűtött szobahőmérséklet melegített

2. faktor: FOOD Táplálékellátottság I I

limitált nem limitált

I

A hernyók tömegét a kikeléstől a bábozódásig mérték.

I

Referencia: J. Kis, F. Kassai, L. Peregovits (nem közölt adatok)

Változók

BOX a dobozok azonosítója, amelyben a hernyókat tartották FOOD Táplálékellátottság limited limitált adlibitum nem limitált TEMPR A fejlődő hernyók környezetének hőmérséklete cooled hűtött room szobahőmérsékletű heated melegített PUPAMASS bábtömeg (g) 1 héttel a bábozódás után STARTMASS kezdeti tömeg (g)

Lepkebáb adatok > pupa[1:5, ] 1 2 3 4 5

BOX 474 423 420 473 430

FOOD limited limited limited limited limited

TEMPR PUPAMASS STARTMASS cooled 0.176 0.007 cooled 0.197 0.013 cooled 0.183 0.015 cooled 0.209 0.015 cooled 0.205 0.017

> str(pupa) ’data.frame’: $ BOX : $ FOOD : $ TEMPR : $ PUPAMASS : $ STARTMASS:

58 obs. of 5 variables: int 474 423 420 473 430 4107 457 80 4101 454 ... Factor w/ 2 levels "adlibitum","limited": 2 2 2 2 2 2 2 2 Factor w/ 3 levels "cooled","heated",..: 1 1 1 1 1 1 1 1 1 num 0.176 0.197 0.183 0.209 0.205 0.197 0.205 0.191 0.234 num 0.007 0.013 0.015 0.015 0.017 0.025 0.034 0.039 0.044

land-use.csv

I

Az adatok egy olyan kutatásból származnak, melyben a mezőgazdaság intenzifikációjának őszi gabona növényekre való hatását vizsgálták a Nagy-Alföldön, szikes talajú területeken.

I

5 gazda, illetve szövetkezet különböző intenzifikációval kezelt földjeit vizsgálták.

I

7 földhasználati intenzitás kategóriát állapítottak meg a felhasznált szerves- és műtrágya valamint növényvédő szer mennyisége alapján (kérdőíves felmérés a vizsgálatok előtt).

land-use.csv

I

3 őszi gabona földet vizsgáltak minden intenzitás szintből.

I

Minden területen kijelöltek 2-2 – 10 db 5 × 1 m-es egymástól 5 m-re egy vonalba eső kvadrátból álló – transzektet, egyet a terület szélén, a másikat pedig 50 m-rel beljebb. Ezzel az elrendezéssel figyelembe vehető a lokális térbeli heterogenitás, ami nagyban meghatározza a biodiverzitási mintázatot.

I

Referencia: Anikó Kovács, Péter Batáry, András Báldi, Andrea Harnos: Weed species richness and cover along an intensification gradient in Central-Hungarian cereal fields (a Weed Research-be beküldve)

Extenzíven használt gabonaföld

Extenzíven használt gabonaföld széle

Intenzíven használt gabonaföld

Intenzíven használt gabonaföld széle

Egy elemzés lépései

I I I I I I I I I

A probléma hátterének megértése. Kérdésfeltevés. Elemző módszerek kiválasztása (adatgyűjtés előtt!) Adatgyűjtés. Az adatok elemezhető formába hozása! Exploratív elemzések (leíró statisztikák, grafikonok). A fő elemzés (modellezés). Interpretáció. Az eredmények közlése.

Milyen egy jó modell?

I I

Tisztáznia kell a dolgokat és nem összezavarni. Parszimóniára kell törekednie. I

I

Things should be made as simple as possible - but no simpler. A. Einstein

Általánosítható. I

az eredményeknek nemcsak a mintánkra kell érvényesnek lennie, hanem arra a statisztikai populációra is, amelyből a megfigyeléseink származnak.

A modellezés folyamata

I I I I I I I I

Kiinduló modell illesztése. A modell redukálása. Modellellenőrzés. Modell megváltoztatása, ha szükséges. Új modell illesztése. Modellredukció. Modellellenőrzés. ...

Az általános lineáris modell General Linear Model

I

I

Egy cél-, vagy függő változó (outcome, response) egyéb változóktól (effects) való függését vizsgáljuk. Az általános lineáris modell speciális esetei: I I I I

Egyszerű Lineáris Regresszió Simple Linear Regression; ANOVA; ANCOVA; stb.

A modell általános felírása I

A lineáris modell: Y = β0 + β1 X1 + β2 X2 + . . . + βp Xp + , Vagy mátrix egyenlet formájában: Y = Xβ + , ahol Y β0 βi X1 , . . . , X p

célváltozó, intercept vagy konstans, ismeretlen paraméterek, magyarázó változók (prediktorok), lehetnek I I

folytonosak (kovariánsok) vagy kvalitatív (faktor) változók. Indikátor (dummy) változók.

A megfigyelt értékeknek és a modell szisztematikus részének a különbsége (mérési hiba vagy nem magyarázott hatás).

A modell praktikus alakja

PUPAMASS = β0 +β1 ·STARTMASS+β2 ·1(TEMPR=high) + + β3 · 1(TEMPR=room) + I

β4 · 1(TEMPR=cooled) hiányzik (túlparaméterezettség)

I

Wilkinson-Rogers formula: PUPAMASS = STARTMASS + TEMPR, ahol a TEMPR egy factor.

Modell formulák R-ben

I I I I

Általános forma: y ∼ tényező1+tényező2+ . . . Interakció a:b, a*b=a+b+a:b Speciális tagok: offset, -1 (nincs konstans) Példák: y ∼ x - Egyszerű Lineáris Regresszió y ∼ x1+x2+x3 - Többszörös Lineáris Regresszió y ∼ f1 - Egytényezős ANOVA y ∼ f1?f2 - Kéttényezős ANOVA interakcióval y ∼ f1?x - ANCOVA

Lepkebáb adatok > pupa[1:5, ] 1 2 3 4 5

BOX 474 423 420 473 430

FOOD limited limited limited limited limited

TEMPR PUPAMASS STARTMASS cooled 0.176 0.007 cooled 0.197 0.013 cooled 0.183 0.015 cooled 0.209 0.015 cooled 0.205 0.017

> str(pupa) ’data.frame’: $ BOX : $ FOOD : $ TEMPR : $ PUPAMASS : $ STARTMASS:

58 obs. of 5 variables: int 474 423 420 473 430 4107 457 80 4101 454 ... Factor w/ 2 levels "adlibitum","limited": 2 2 2 2 2 2 2 2 Factor w/ 3 levels "cooled","heated",..: 1 1 1 1 1 1 1 1 1 num 0.176 0.197 0.183 0.209 0.205 0.197 0.205 0.191 0.234 num 0.007 0.013 0.015 0.015 0.017 0.025 0.034 0.039 0.044

Változók BOX a dobozok azonosítója, amelyben a hernyókat tartották FOOD Táplálékellátottság limited limitált adlibitum nem limitált TEMPR A fejlődő hernyók környezetének hőmérséklete cooled hűtött room szobahőmérsékletű heated melegített PUPAMASS bábtömeg (g) 1 héttel a bábozódás után STARTMASS kezdeti tömeg (g)

Exploratív elemzések 1. lépés: Adatok leíró statisztikái és grafikonjai. Numerikus áttekintés: > summary(pupa) BOX Min. : 39.0 1st Qu.: 409.2 Median : 443.0 Mean :1015.5 3rd Qu.: 473.8 Max. :4202.0 STARTMASS Min. :0.00100 1st Qu.:0.00700 Median :0.01400 Mean :0.01895 3rd Qu.:0.02375 Max. :0.11900

FOOD adlibitum:32 limited :26

TEMPR cooled:19 heated:22 room :17

PUPAMASS Min. :0.1430 1st Qu.:0.1970 Median :0.2560 Mean :0.2571 3rd Qu.:0.3103 Max. :0.4330

Grafikus sűrűségfüggvény becslések

I

I

Hisztogram: durva becslés - érzékeny az osztályintervallumok megválasztására. Kernel sűrűség becslés: jobb. Simított hisztogram.

Hisztogram >hist(PUPAMASS,main="Bábtömeg egy héttel a bábozódás után",xlab="Tömeg (g)") Bábtömeg egy héttel a bábozódás után

Frequency

15

10

5

0 0.10

0.15

0.20

0.25

0.30

Tömeg (g)

0.35

0.40

0.45

Simított hisztogram >plot(density(PUPAMASS),main="Bábtömeg egy héttel a bábozódás után"); rug(PUPAMASS)

3 2 1 0

Density

4

5

Bábtömeg egy héttel a bábozódás után

0.1

0.2

0.3 N = 58 Bandwidth = 0.02654

0.4

0.5

Boxplotok

I I I

Kvalitatív és kvantitatv változók kapcsolata. Faktor kombinációk is lehetnek a vízszintes tengelyen. Különbségek és interakciók becslése.

Boxplot >boxplot(PUPAMASS ∼(TEMPR:FOOD),col=2:4, names=c("CA","HA","RA","CL","HL","RL"))

0.40

●

0.25

0.30

0.35

● ●

●

● ●

0.15

0.20

●

● ●

CA

HA

RA

CL

HL

RL

Hegedűábra

●

0.20

0.25

0.30

0.35

0.40

vioplot(PUPAMASS[FOOD=="adlibitum"], PUPAMASS[FOOD=="limited"], col="white", names=c("ad libitum", "limited"))

0.15

●

ad libitum

limited

Interakciós ábra

I

I

A célváltozó különböző faktor kombinációk szerinti átlagait ill. más leíró statisztikáit rajzolja ki így ábrázolva a lehetséges interakciókat (nem additív hatásokat). Ha a vonalak többé-kevésbé párhuzamosak, akkor nem várunk interakciót.

Interakciós ábra

FOOD

0.20 0.22 0.24 0.26 0.28 0.30

mean of PUPAMASS

>interaction.plot(TEMPR,FOOD,PUPAMASS)

adlibitum limited

cooled

heated

room TEMPR

Szórásdiagram >plot(PUPAMASS ∼ STARTMASS,main="PUPAMASS-STARTMASS Szórásdiagram", xlab="STARTMASS",ylab="PUPAMASS",pch=".") PUPAMASS−STARTMASS Szórásdiagram

0.40

●

0.35

● ● ●

● ●

0.25

0.30

●

● ● ● ● ● ●●● ● ●● ● ● ● ●

●

● ●

●

● ●

●

●

● ●

0.20

● ●

●

●

●

●

● ●

●

●

● ●

●

●

●

●

●

● ●

●

●

0.15

PUPAMASS

●

●

●

●

●

0.00

0.02

0.04

0.06 STARTMASS

0.08

0.10

0.12

Feltételes szórásdiagram >coplot(PUPAMASS STARTMASS|TEMPR*FOOD, xlab="STARTMASS",ylab="PUPAMASS",pch=20) Given : TEMPR room heated cooled 0.08

0.12

0.00

0.04

0.08

0.12

● ●

●

● ●

● ● ● ●

● ●

● ●

● ●

●

●

0.35

●

●

●

● ● ● ● ● ● ●

●

●

● ● ●

● ●

●

● ●

●

●

●

0.00

0.04

0.08

STARTMASS

0.12

adlibitum

●

● ● ● ● ● ●

0.25

limited

0.25

●

●

0.15

●

0.15

PUPAMASS

● ●

● ●● ● ●

●

Given : FOOD

0.04

0.35

0.00

Egyszerű Lineáris Regresszió 1. modell: PUPAMASS = β0 + β1 STARTMASS + > mod1.lm <- lm(PUPAMASS ~ STARTMASS) > summary(mod1.lm) Call: lm(formula = PUPAMASS ~ STARTMASS) Residuals: Min 1Q Median 3Q -0.115249 -0.053640 -0.003456 0.047954

Max 0.169238

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.26866 0.01217 22.071 <2e-16 STARTMASS -0.61265 0.45143 -1.357 0.180 Residual standard error: 0.06594 on 56 degrees of freedom Multiple R-Squared: 0.03184, Adjusted R-squared: 0.01455 F-statistic: 1.842 on 1 and 56 DF, p-value: 0.1802

Szórásdiagram az illesztett egyenessel >abline(mod1.lm) PUPAMASS−STARTMASS szórásdiagram

0.40

●

0.35

● ● ●

● ●

0.25

0.30

●

● ● ● ● ● ●●● ● ●● ● ● ● ●

●

● ●

●

● ●

●

●

● ●

0.20

● ●

●

●

●

●

● ●

●

●

● ●

●

●

●

●

●

● ●

●

●

0.15

PUPAMASS

●

●

●

●

●

0.00

0.02

0.04

0.06 STARTMASS

0.08

0.10

0.12

Kéttényezős ANOVA

> mod2.lm <- lm(PUPAMASS ~ FOOD + TEMPR) > anova(mod2.lm) Analysis of Variance Table Response: PUPAMASS Df Sum Sq Mean Sq F value Pr(>F) FOOD 1 0.161984 0.161984 98.906 8.306e-14 TEMPR 2 0.001065 0.000532 0.325 0.724 Residuals 54 0.088439 0.001638

Kéttényezős ANOVA interakcióval

> mod3.lm <- lm(PUPAMASS ~ FOOD * TEMPR) > anova(mod3.lm) Analysis of Variance Table Response: PUPAMASS Df Sum Sq Mean Sq F value Pr(>F) FOOD 1 0.161984 0.161984 99.7353 1.088e-13 TEMPR 2 0.001065 0.000532 0.3277 0.7220 FOOD:TEMPR 2 0.003984 0.001992 1.2264 0.3017 Residuals 52 0.084455 0.001624

Kéttényezős ANOVA interakcióval és kovariánssal

> mod4.lm <- lm(PUPAMASS ~ FOOD * TEMPR + STARTMASS) > anova(mod4.lm) Analysis of Variance Table Response: PUPAMASS Df Sum Sq Mean Sq F value Pr(>F) FOOD 1 0.161984 0.161984 97.8787 1.857e-13 TEMPR 2 0.001065 0.000532 0.3216 0.7264 STARTMASS 1 0.000298 0.000298 0.1804 0.6729 FOOD:TEMPR 2 0.003738 0.001869 1.1294 0.3312 Residuals 51 0.084402 0.001655

Hipotézisvizsgálatok

I

I

A modell egy vagy több prediktorának szignifikanciáját állapítjuk meg. Ha a hibatagok I I

I

függetlenek és normális eloszlásúak.

Két beágyazott modell (a szűkebb modell magyarázó változóinak halmaza részhalmaza a bővebb modellének) összehasonlítható egy F-teszttel: anova(model1, model2).

Hipotézisvizsgálatok

I

I

I

Egy általánosan végzett tesz az aktuális modell nullmodellhez való hasonlítása (nincsenek prediktorok, csak a konstans (intercept)): anova(model). (A modell egészének szignifikanciája.) Egyedi prediktorok F-próbával, vagy egy t-próbával summary(model) tesztelhetők. Kerüljük a t-tesztek használatát kettőnél több szintű kvalitatív predictorok (faktorok) esetén!

A modellezés céljai

I

Predikció: I

I

I

I

A változók közötti kapcsolat megértése. I

I I

Megfigyelünk új X-eket és a hozzá tartozó Y -t szeretnénk megbecsülni. A predikciós teljesítmény javul az olyan változók eltávolításáaval, amik nem nagyon játszanak szerepet. Automatikus változó szelekciók jól működhetnek. Manuális szelekció jobb.

Gyakran mindkettő célunk. Nam tanácsos teljesen automatikus szelekciós módszerekre hagyatkozni.

Változó szelekciós módszerek I.

I

Akaike Information Criterion (AIC) AIC = −2logLik + 2p,

I I

ahol p a paraméterek száma. Általános, normál lineáris modelleken túl is használható. A step() függvény ezt használja. I I I

Lépésenkénti keresés a lehetséges modellek terében. Szekvenciálisan távolít el (vagy vesz be) változókat. Minimalizálja az AIC-ot.

Változó szelekciós módszerek II.

I

Tesztelésen alapulnak. F -teszttel hasonlítják össze a beágyazott modelleket. Nem igazán jó módszer: a beválasztott változók sorrendje nagyon számít. Rosszabb, mint a kritériumra épülő módszerek. Manuális változó szelekcióra használható.

I

drop1(model,test="F")

I I I

I

Automatikus változó szelekció

> mod5.lm <- step(mod4.lm, trace = 0) > anova(mod5.lm) Analysis of Variance Table Response: PUPAMASS Df Sum Sq Mean Sq F value Pr(>F) FOOD 1 0.161984 0.161984 101.35 3.585e-14 Residuals 56 0.089503 0.001598

Két modell összehasonlítása

> anova(mod4.lm, mod5.lm) Analysis of Variance Table Model 1: PUPAMASS ~ FOOD * TEMPR + STARTMASS Model 2: PUPAMASS ~ FOOD Res.Df RSS Df Sum of Sq F Pr(>F) 1 51 0.084402 2 56 0.089503 -5 -0.005101 0.6165 0.6877 I I

Nincs szignifikáns különbség a modellek között. Válasszuk a szűkebb modellt!

Manuális változó szelekció

> drop1(mod4.lm, test = "F") Single term deletions Model: PUPAMASS ~ FOOD * TEMPR + STARTMASS Df Sum of Sq RSS AIC F value Pr(F) <none> 0.08 -364.89 STARTMASS 1 5.291e-05 0.08 -366.85 0.0320 0.8588 FOOD:TEMPR 2 3.738e-03 0.09 -366.38 1.1294 0.3312

Konfidencia-intervallumok

I I I

Tartományok a paraméterek lehetséges értékeire. A hatásnagyságok becslésére hasznosabb, mint a p-érték. A p-értékek a statisztikai szignifikanciát mutatják, nem pedig a gyakorlati jelentőséget.

> confint(mod4.lm) (Intercept) FOODlimited TEMPRheated TEMPRroom STARTMASS FOODlimited:TEMPRheated FOODlimited:TEMPRroom

2.5 % 97.5 % 0.28712994 0.33963231 -0.15362149 -0.07630159 -0.04029847 0.03192612 -0.06188462 0.01129474 -0.54923885 0.65664283 -0.05702646 0.04739676 -0.02174418 0.08834948

Diagnosztika

A lineáris modell feltételeinek ellenőrzése. I Korrekt-e a modell szisztematikus része (linearitás)? I A modell véletlen részét () tekintve: I I I

I

konstans variancia, korrelálatlanság, normalitás.

Torzító pontok keresése (olyan pontok, amelyeknek a többi pontnál sokkal nagyobb hatása van az illesztett modellre).

Diagnosztikus módszerek I I

Lehetnek numerikusak vagy grafikusak. Általában a grafikus módszereket preferáljuk, mert informatívabbak. I I

I

I

I

I

reziduális ábrák, normalitást ellenőrző ábrák.

Gyakorlatilag lehetetlen megállapítani egy modellről, hogy teljesen korrekt-e. A diagnosztikák célja: leellenőrizni, hogy a modell nem durván rossz-e. Több figyelmet kell fordítani arra, hogy ne kövessünk el nagy hibákat, mint arra, hogy a modellünk optimális-e. Négy hasznos ábra: plot(model)

● ●● ● ●

●

● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●

● ● ●

●●

● ●

0.20

0.22

0.24

0.26

43

0.28

0.30

0.32

Normal Q−Q 33 ● 39 ●

2

● ●

0

0.00

●● ● ● ●● ●

● ● ● ● ● ● ●● ●

39

● ● ● ● ●

● ● ●●

●

−2

● ●

0.22

0.24

0.26

Fitted values

0.28

● ●

−1

0

1

2

● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●

0.30

0.32

Residuals vs Leverage ●

33

39 ●

2

33 ●

● ● ● ● ● ●● ● ● ●● ●● ●● ●● ● ●● ● ●●● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ●● ● ● ●● ● ● ●

0

● ● ● ● ● ● ● ●

● ●

0.20

43 39

1 0.5 ● ●

0.5 1

●

Cook's distance ●

−4

● ● ●

● ●

Standardized residuals

1.0 0.0


● ● ●● ●●

●

Theoretical Quantiles

Scale−Location

● ● ● ●● ●● ● ●

● ●

43

Fitted values

● ● ●

●● ●● ●●●● ●●●● ●●●●●●●●●● ●●●●●●●● ●●●●●●● ●●●●●● ● ●●●

●

−3

33 ● ●

−0.15

Residuals

0.15

Residuals vs Fitted ●●


Diagnosztikus ábrák

43

0.0

0.1

0.2

0.3 Leverage

0.4

0.5

0.6

Illeszkedés ellenőrzése >plot(mod4.lm,1,pch=20)

0.15

Residuals vs Fitted

●

● ●

● ●● ●● ● ●

−0.05

39 ●

● ●

● ●● ●● ●

● ● ● ●

● ●

●

●

● ● ● ●

● ● ● ●● ● ●

● ● ●

● ●

● ● ● ●● ● ● ● ●

● ●

●

−0.15

Residuals

0.05

0.10

33 ● ●

●

0.20

0.22

0.24

0.26

43

0.28

Fitted values lm(PUPAMASS ~ FOOD * TEMPR + STARTMASS)

0.30

0.32

Reziduumok normalitása >plot(mod4.lm,2,pch=20) Normal Q−Q 3

33 ●

2 0

1

●

−1

● ● ● ●

−2

●

●

●

●●● ● ●

●●●

●●●

●●●

● ●●●●

●●●● ●

●●● ●

●●● ●

●●●

● ●● ●

●

●

●

●

●

●

−3


39 ● ● ●

●

43

−2

−1

0

1

Theoretical Quantiles lm(PUPAMASS ~ FOOD * TEMPR + STARTMASS)

2

A variancia állandóságának ellenőrzése >plot(mod4.lm,3,pch=20) Scale−Location ●

43

1.5

●

39 ● ●

●

1.0

● ●

●

● ●

● ●

● ●

●

●

● ●

0.5

● ● ● ●

● ●

● ●

●

● ●

●● ● ●

● ● ●

●● ●

● ● ● ● ●

●

● ●

●

●

●

● ●

● ● ● ●

0.0


33 ●

0.20

0.22

0.24

0.26

0.28

Fitted values lm(PUPAMASS ~ FOOD * TEMPR + STARTMASS)

0.30

0.32

Torzító pontok keresése >plot(mod4.lm,5,pch=20) Residuals vs Leverage 33

●

2

1 ●

● ●

●●

0.5 ●

●

●● ● ● ●

●

0

● ●● ● ● ●● ●● ● ●● ● ● ● ●● ● ●● ● ● ● ● ●● ● ●● ● ● ●

−2

● ●

●

● ●

● ●

●

●

0.5

●

1

●

43

Cook's distance

−4


39 ●

0.0

0.1

0.2

0.3

0.4

Leverage lm(PUPAMASS ~ FOOD * TEMPR + STARTMASS)

0.5

0.6

Hogy detektáljuk a problémákat?

I

I

I

Illeszkedés ellenőrzése: A reziduumokban nem lehet trend (y = 0). Ha van, meg kell változtatni a modellt (transzformáció, nemlineáris modell etc). Reziduumok normalitása: QQ-ábra. A reziduumokat az "ideális" normális eloszlású megfigyelésekhez hasonlítjuk. Normális eloszlás esetén a pontok lineáris trendet követnek (y = x). Egyébként ferdeséget jeleznek. Scale-location ábra: a variancia homogenitását lehet vele ellenőrizni.

Hogy detektáljuk a problémákat?

I

I

Residuals vs. Leverage ábra: Torzító pontok keresése. A pontoknak az adott Cook távolság (Cook’s distance) szinteken belül kell lennie. A számozott pontok lehetnek gyanúsak. Cook-féle távolság: az illeszkedés megváltozásának standardizált mértéke, ha az adott megfigyelést kivesszük az adatok közül.

Cook’s distance plot >plot(mod4.lm,4,pch=20)

0.15

43

0.10

33

0.05

39

0.00

Cook's distance

0.20

Cook's distance

0

10

20

30

40

Obs. number lm(PUPAMASS ~ FOOD * TEMPR + STARTMASS)

50

60

Lineáris modell - korlátok

Nagyon sok kapcsolatot nem írható le egyszerű lineáris modellel, mivel I a függő változó lehet nem folytonos (és nem normális) eloszlású (pl. gyakoriságok, bináris adatok); I a magyarázó változók hatása a függő változóra lehet, hogy nem lineáris; I a megfigyelési egységek lehet, hogy nem függetlenek; I a variancia lehet, hogy nem konstans.

Általánosított lineáris modellek (Generalized Linear Models)

Az általános lineáris modell általánosítása: I Megengedi, hogy az eloszlás nem normális legyen (pl. Poisson, binomiális ill., multinomiális (exponenciális eloszláscsalád)). I A variancia állandóságának feltétele sem olyan szigorú, mint a hagyományos lineáris modelleknél.

Hogy általánosít ez a módszer? I

I

A függő változót most is a magyarázó változók lineáris kombinációjából becsüljük. A függő és magyarázó változók egy ún. „link” függvénnyel vannak összekapcsolva: η = β0 + β1 X1 + β2 X2 + . . . + βk Xk , lineáris egyenlet, ahol η lineáris prediktor, X magyarázó változók, β együtthatók.

I

Maximum likelihood (ML) módszerrel illesztünk. g(Y ) = η link függvény.

I

glm(formula, family = gaussian, ...)

I

Súgó a függvény családról: ?family

I

Gyakorisági adatok regressziója (count regression)

I I

I

I

I

A függő változó gyakorisági adat (pozitív egész). Ha az összes lehetőség egy adott korlátos szám, akkor binomiális modellt használunk. Van-nincs (0-1) adatok esetén a binomiális modell használatos (logisztikus regresszió). Ha a gyakoriságok elegendően nagyok, akkor az általános lineáris modell is jó lehet. Egyéb esetekben a Poisson és - kevésbé gyakran - a negatív binomiális modell használható.

Poisson regresszió

Ha Y Poisson eloszlású µ > 0 várható értékkel, akkor: P (Y = y) = E(Y ) = var(Y ) = µ.

eµ µy , y!

y = 0, 1, 2, ...

Honnan származhatnak Poisson-eloszlású adatok? I

I

I

Ha a gyakoriságok egy előre rögzített számú megfigyelésből származnak, akkor a függő változót binomiálisként modellezhetjük. Kis siker valószínűségek, és nagyszámú összes lehetőség esetén alkalmazhatjuk a Poisson közelítést. (Pl. ritka incidenciája egy adott fajnak egy földrajzi területen.) Ha gyakoriságokat számolunk egy adott időintervallumban, területen, térrészben, anyagmennyiségben, és a siker valószínűsége arányos az intervallum hosszával, térrész térfogatával stb., és független más eseményektől. (Pl. bejövő telefonhívások, földrengések száma stb.) Fontos: Poisson-eloszlású véletlen változók összege is Poisson. (Hasznos, ha csak aggregált adataink vannak.)

Földhasználati példa 36 mintavételi terület esetén vannak adataink a következőkről: Weedcover Az adott transzekt teljes gyomborítottsága százalékosan. Totspeciesnb A gyomnövény fajok száma. N input Éves nitrogén bevitel. Transectpos A transzekt elhelyezkedése a földterületen. I 0 - a transzekt közvetlenül a terület szélén helyezkedik el, I 1 - belül van. Transect pair Ugyanahhoz a földhöz tartozó transzekt párok azonosítója.

Földhasználati példa

Noncrop area A tanulmányozott transzekt körül húzott 500 m sugarú körbe eső nem művelt terület százalékos aránya (főleg füves terület, de lehet erdős, beépített, mocsaras vagy nyílt vizes terület). Modellezni szeretnénk a gyomnövény fajok számát és a gyomborítottságot a nitrogén bevitel, a nem művelt terület aránya és a transzekt pozíció függvényében.

Az adatok struktúrája

> str(land) ’data.frame’: 42 obs. of 10 variables: $ SampleArea : Factor w/ 42 levels "AG30E ","AG30I ",..: 31 32 33 34 $ Weedcover : int 53 44 48 29 61 45 33 33 36 38 ... $ Totspeciesnb : int 29 21 37 28 29 21 36 29 38 13 ... $ Intensity : int 1 1 1 1 1 1 2 2 2 2 ... $ N_input : int 0 0 0 0 0 0 113 113 113 113 ... $ Herbicide_use: int 0 0 0 0 0 0 1 1 1 1 ... $ Transectpos : Factor w/ 2 levels "0","1": 1 2 1 2 1 2 1 2 1 2 ... $ Transect_pair: int 16 16 17 17 18 18 19 19 20 20 ... $ Noncrop_area : int 49 50 41 43 55 53 21 21 9 8 ... $ Farmer : Factor w/ 5 levels "AG","ET","NL",..: 5 5 5 5 5 5 5 5

Density plot >plot(density(Totspeciesnb))

0.02 0.01 0.00

Density

0.03

density.default(x = Totspeciesnb)

0

10

20

30 N = 42 Bandwidth = 4.453

40

50

60

Boxplot

30 20 10

Totspeciesnb

40

50

>plot(Totspeciesnb ∼ Transectpos)

0

1 Transectpos

Hegedűábra

●

20

30

40

50

>vioplot(Totspeciesnb[Transectpos==0],Totspeciesnb[Transectpos== col="white")

10

●

1

2

Szórásdiagram >plot(Totspeciesnb ∼ Noncrop area)

50

●

●

●

40

●

● ●

●

●

30

●

●

●

●

●

● ●

●

● ●

●

20

●

●

●

● ● ● ● ● ●

●

●

●

●

● ●

● ●

●

●

10

Totspeciesnb

●

●

●

●

0

10

20

30 Noncrop_area

40

50

60

Feltételes szórásdiagram >coplot(Totspeciesnb Noncrop area|Transectpos,pch=20) Given : Transectpos

1

0

0

10

20

30

40

50

60

50

●

40 30

●

●

●

●

●

●

● ●

● ●

●

●

●

●

●

●

●

20

●

●

●

● ●

●

● ●

●

●

●

● ● ●

10

Totspeciesnb

●

●

●

●

●

●

●

●

●

●

0

Noncrop_area

10

20

30

40

50

60

Interakciós ábra

45

>interaction.plot(as.factor(N-input),Transectpos,Totspeciesnb)

40

Transectpos

30 25 20 15

mean of Totspeciesnb

35

0 1

0

34

68

92

100

as.factor(N_input)

113

270

Lineáris modell > mod1.lm <- lm(Totspeciesnb ~ N_input * Transectpos + + Noncrop_area) > summary(mod1.lm) Call: lm(formula = Totspeciesnb ~ N_input * Transectpos + Noncrop_area) Residuals: Min 1Q Median 3Q Max -13.5241 -4.7578 -0.3224 4.1886 17.9341 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 29.88544 3.68671 8.106 1.01e-09 N_input -0.03013 0.02115 -1.424 0.162731 Transectpos1 -15.68275 3.60955 -4.345 0.000104 Noncrop_area 0.13763 0.06801 2.024 0.050266 N_input:Transectpos1 0.01105 0.02881 0.384 0.703520 Residual standard error: 7.437 on 37 degrees of freedom Multiple R-Squared: 0.5816, Adjusted R-squared: 0.5364 F-statistic: 12.86 on 4 and 37 DF, p-value: 1.172e-06

Reziduum vs. becsült érték ábra >plot(mod1.lm,1,pch=20)

20

Residuals vs Fitted

33 ●

8

●

10

●

●

●

●

● ●

●

●

● ●

●

●

0

●

● ●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

● ●

●

●

−10

Residuals

● ● ●

● ●

●

31 ●

10

15

20

25

30

Fitted values lm(Totspeciesnb ~ N_input * Transectpos + Noncrop_area)

35

Normalitás vizsgálat (Normal QQ-plot) >plot(mod1.lm,2)

3

Normal Q−Q

2

33 ●

●

8●

●

1

●

●

● ● ● ● ●

●

● ● ● ●

0

● ● ●

● ●

● ● ● ●

● ●

● ●

●

●

−1

●

●

●

●

●

●

●

● ●

−2


●

●

31

−2

−1

0

1

Theoretical Quantiles lm(Totspeciesnb ~ N_input * Transectpos + Noncrop_area)

2

Szórás-becsült érték ábra (Scale-location plot) >plot(mod1.lm,3,pch=20) Scale−Location 1.5

33 ●

●

8

31 ●

●

● ●

●

● ● ●

1.0

●

●

● ●

● ●

●

●

● ●

●

● ● ●

●

●

● ●

●

● ●

0.5

●

● ● ● ●

● ●

●

● ●

0.0


●

10

15

20

25

30

Fitted values lm(Totspeciesnb ~ N_input * Transectpos + Noncrop_area)

35

Problémák

I I I I

Enyhén nemlineáris trend. Nem konstans variancis. Enyhén nem normális eloszlású hibatag. Próbáljuk meg transzformálni az adatokat, pl. logaritmus transzformáció!

Lineáris modell log transzformált függő változóval > mod2.lm <- lm(log(Totspeciesnb + 1) ~ N_input * Transectpos + + Noncrop_area) > summary(mod2.lm) Call: lm(formula = log(Totspeciesnb + 1) ~ N_input * Transectpos + Noncrop_area) Residuals: Min 1Q Median 3Q Max -0.55314 -0.20587 0.02715 0.18937 0.64625 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.406e+00 1.422e-01 23.945 < 2e-16 N_input -1.374e-03 8.161e-04 -1.684 0.100570 Transectpos1 -6.005e-01 1.393e-01 -4.312 0.000115 Noncrop_area 4.777e-03 2.624e-03 1.821 0.076760 N_input:Transectpos1 4.096e-05 1.111e-03 0.037 0.970800 Residual standard error: 0.2869 on 37 degrees of freedom Multiple R-Squared: 0.6147, Adjusted R-squared: 0.5731 F-statistic: 14.76 on 4 and 37 DF, p-value: 2.685e-07

Reziduum vs. becsült érték ábra >plot(mod2.lm,1,pch=20) Residuals vs Fitted

0.6

●

8

●

●

●

●

0.2

●

● ●

●

●

● ●

●

0.0

● ●

●

●

●

●

●

● ● ●●

●

−0.2

● ●

●

● ●

●

● ● ●

●

−0.4

● ● ●

●

31 ●

−0.6

Residuals

0.4

33 ● ●

2.6

2.8

3.0

3.2

3.4

Fitted values lm(log(Totspeciesnb + 1) ~ N_input * Transectpos + Noncrop_area)

3.6

Normal QQ-plot >plot(mod2.lm,2,pch=20) Normal Q−Q

2

8●

42 ● ●

1

● ●

●

●

●

●

●

●

●

● ● ●

0

●

● ● ● ● ●

● ● ●

●

−1

● ●

●

●

●

● ●

● ●

●

●

● ● ● ●

−2


●

●

31

−2

−1

0

1

Theoretical Quantiles lm(log(Totspeciesnb + 1) ~ N_input * Transectpos + Noncrop_area)

2

Scale-location plot >plot(mod2.lm,3,pch=20)

1.5

Scale−Location ●

8 31 ●

●

42 ●

●

● ● ●

●

●

1.0

● ●

●

●

●

●

●

●

●

●

● ●

● ●

●

●

●

●

● ● ●

●

0.5

●

●

●

●

●

●

● ●

●

0.0


●

2.6

2.8

3.0

3.2

3.4

Fitted values lm(log(Totspeciesnb + 1) ~ N_input * Transectpos + Noncrop_area)

3.6

A két modell összehasonlítása

> summary(mod1.lm)$adj.r.squared [1] 0.5364204

> summary(mod2.lm)$adj.r.squared [1] 0.573088 I I I

Nem nagy javulás. Jobb illeszkedés. Nehézkes interpretáció.

Poisson modell

> mod1.pois <- glm(Totspeciesnb ~ N_input * Transectpos + + Noncrop_area, family = poisson) > mod1.pois

Call: glm(formula = Totspeciesnb ~ N_input * Transectpos + Noncrop_area Coefficients: (Intercept) N_input Transectpos1 3.3428928 -0.0009255 -0.5800613 Noncrop_area N_input:Transectpos1 0.0057278 -0.0006015 Degrees of Freedom: 41 Total (i.e. Null); Null Deviance: 198.2 Residual Deviance: 76.93 AIC: 293.5

37 Residual

summary(mod1.pois) Call: glm(formula = Totspeciesnb ~ N_input * Transectpos + Noncrop_area, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -2.6342 -1.1127 -0.1734 0.8111 3.1884 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.3428928 0.0975010 34.286 < 2e-16 N_input -0.0009255 0.0005436 -1.702 0.08868 Transectpos1 -0.5800613 0.1018254 -5.697 1.22e-08 Noncrop_area 0.0057278 0.0018755 3.054 0.00226 N_input:Transectpos1 -0.0006015 0.0008987 -0.669 0.50329 (Dispersion parameter for poisson family taken to be 1) Null deviance: 198.232 Residual deviance: 76.934 AIC: 293.51

on 41 on 37

degrees of freedom degrees of freedom

Number of Fisher Scoring iterations: 4

Hatások tesztelése

I

I I

I

A summary(model) közelítő Wald teszteket csinál. Az SE-k lehet, hogy túlbecsültek, és így elvesztünk szignifikáns eredményeket. A deviancia alapú tesztek jobbak. A deviancia azt méri, hogy a modell mennyire van közel a tökéleteshez. (A lineáris modell esetén: deviancia = RSS.) Chi2 eloszlású. A determinációs együttható (R-négyzet a lineáris modelleknél): > 1 - 77/198 [1] 0.6111111

Anova a Poisson modellre > anova(mod1.pois, test = "Chi") Analysis of Deviance Table Model: poisson, link: log Response: Totspeciesnb Terms added sequentially (first to last)

NULL N_input Transectpos Noncrop_area N_input:Transectpos

Df Deviance Resid. Df Resid. Dev P(>|Chi|) 41 198.232 1 15.488 40 182.744 8.302e-05 1 96.027 39 86.717 1.133e-22 1 9.332 38 77.384 0.002 1 0.451 37 76.934 0.502

Poisson modell interakciókkal > mod2.pois <- glm(Totspeciesnb ~ (N_input + Noncrop_area + + Transectpos)^2, family = poisson) > anova(mod2.pois, test = "Chi") Analysis of Deviance Table Model: poisson, link: log Response: Totspeciesnb Terms added sequentially (first to last)

NULL N_input Noncrop_area Transectpos N_input:Noncrop_area N_input:Transectpos Noncrop_area:Transectpos

Df Deviance Resid. Df Resid. Dev P(>|Chi|) 41 198.232 1 15.488 40 182.744 8.302e-05 1 9.914 39 172.829 0.002 1 95.445 38 77.384 1.521e-22 1 1.016 37 76.368 0.313 1 0.462 36 75.906 0.497 1 3.123 35 72.783 0.077

A két modell összehasonlítása

> anova(mod1.pois, mod2.pois, test = "Chi") Analysis Model 1: Model 2: Resid. 1 2

of Deviance Table Totspeciesnb ~ N_input * Transectpos + Noncrop_area Totspeciesnb ~ (N_input + Noncrop_area + Transectpos)^2 Df Resid. Dev Df Deviance P(>|Chi|) 37 76.934 35 72.783 2 4.151 0.126

Modell szelekció

> drop1(mod2.pois, test = "Chi") Single term deletions Model: Totspeciesnb ~ (N_input + Noncrop_area + Transectpos)^2 Df Deviance AIC LRT Pr(Chi) <none> 72.783 293.359 N_input:Noncrop_area 1 73.747 292.323 0.964 0.32607 N_input:Transectpos 1 74.329 292.905 1.546 0.21367 Noncrop_area:Transectpos 1 75.906 294.482 3.123 0.07719

Diagnosztikus ábrák

1

● ●

●

●

● ●

● ●

●

● ●

●

●

●

●

●

●

●●

●●

●

●

●

● ● ●

●

●

● ●

●

●

●

●

3 2

●●●

−3

−3

2.6

2.8

3.0

3.2

3.4

●

●

31 ●

2.4

3.6

●

31

● ● ● ●

●

−2

●

●

●●

●

● ●●

0

1

2

4 ●

●

●

●

●

●●

●

● ●

●

●

●

●

●

●

●●

●

● ●

●

●

●

●

●

● ●

2.8

3.0

3.2

Predicted values

3.4

3.6

●

●

0.5

42

3

● ● ●

● ●

● ●

●

●

● ●

●

●

●

● ●● ● ● ●

●

● ●

●

● ● ●

●

● ●

● ●

41 ●

●

●

0.5 39 ●

● Cook's distance

−3

●

1 ●

2

●

●

●

1.0

●

●

1

● ●

0.5

●●●

●●●

●

−1

●

●

Std. deviance resid.

1.5

33 ●

31 ●

0.0


●

●●

●

Residuals vs Leverage

8 ●

2.6

●

●●●

−1

Scale−Location

2.4

●

Theoretical Quantiles

●

●

●

●

●

Predicted values

●

8●

33 ● ● ●

1

●

● ●

−1 0

3 2

33 ●

●

−1 0

Residuals

Normal Q−Q

8


Residuals vs Fitted ●

0.00

0.05

0.10

0.15

0.20

Leverage

0.25

1

0.30

Illeszkedés ellenőrzése >plot(mod1.pois,1,pch=20) Residuals vs Fitted

8

3

●

33 ●

●

2

● ● ● ●

1

●

●

● ● ●

●

0

●

●

●

●

●

●

●●

●

−1

●

●

●

● ● ●

●

●

●

●

● ●

●

−2

●

●

31 ●

−3

Residuals

●

● ● ●

2.4

2.6

2.8

3.0

3.2

Predicted values glm(Totspeciesnb ~ N_input * Transectpos + Noncrop_area)

3.4

3.6

Parciális reziduális ábrák >library(gam) > par(mfrow=c(1,3),pty="s") >plot.gam(mod1.pois,resid=T,pch=20) >par(mfrow=c(1,1))

0

1

3

●

3

●

3

●

●

●

●

●●

● ●

●

● ●

● ●

● ●

−1

●

● ●

●

●

●

●

50

100

2 ● ● ● ●

● ● ●● ● ●

● ● ● ● ●●

●

150 N_input

200

250

1

●

● ●

●

●

●

● ● ● ● ●

● ●

●

● ●

● ●

●

●● ●

●

●

●

●

Transectpos

●

● ● ●

● ●●

● ●

● ●

● ●

●

●

●

0

● ●

●

●

●

●

●

●

●

●

0

●

● ●

●

●

−1

●

●

●

●

partial for Noncrop_area

0

●

● ●

●

●

−2

● ●

1

●

−2

1

●

●

● ●●

0

●

●

●

partial for Transectpos

● ●

−2

partial for N_input

●

●

●

2

●

−1

2

●

●

0

10

20

30

40

Noncrop_area

50

60

Illeszkedés ellenőrzése

I

I

Az ábrák majdnem ugyanúgy használhatók, mint a lineáris modell esetén. A normalitás általában nem teljesül tökéletesen. A parciális reziduálisok ellenőrzésére a plot.gam használható a gam csomagból.

Túlszóródás I I

I I

I

I

Poisson változó esetén az átlag és a variancia megegyezik. A variancia függvényt az átlag teljesen meghatározza, nem szabad paraméter. Az ún. diszperziós paraméter 1. Gyakran túlságosan szigorú ez a feltétel. Gyakran túlszóródás (overdispersion) van. A túlszóródást a reziduális deviancia és a hozzá tartozó szabadsági fokból határozható meg. Többé-kevésbé egynelőnek kell lenniük. Ha nagyon különbözőek, akkor az ún. quasilikelihood módszert használhatjuk, amellyel a modellparaméterek a hiba eloszlás teljes ismerete nélkül határozhatók meg.

A diszperziós paraméter ellenőrzése

> deviance(mod1.pois)/df.residual(mod1.pois) [1] 2.079286

Poisson modell túlszóródással > mod1.qpois <- glm(Totspeciesnb ~ N_input * Transectpos + + Noncrop_area, family = quasipoisson) > summary(mod1.qpois) Call: glm(formula = Totspeciesnb ~ N_input * Transectpos + Noncrop_area, family = quasipoisson) Deviance Residuals: Min 1Q Median 3Q Max -2.6342 -1.1127 -0.1734 0.8111 3.1884 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.3428928 0.1426271 23.438 < 2e-16 N_input -0.0009255 0.0007952 -1.164 0.251962 Transectpos1 -0.5800613 0.1489530 -3.894 0.000398 Noncrop_area 0.0057278 0.0027436 2.088 0.043767 N_input:Transectpos1 -0.0006015 0.0013147 -0.458 0.649948 (Dispersion parameter for quasipoisson family taken to be 2.139865) Null deviance: 198.232 Residual deviance: 76.934

on 41 on 37

degrees of freedom degrees of freedom

Poisson modell túlszóródással > anova(mod1.qpois, test = "F") Analysis of Deviance Table Model: quasipoisson, link: log Response: Totspeciesnb Terms added sequentially (first to last)

NULL N_input Transectpos Noncrop_area N_input:Transectpos

Df Deviance Resid. Df Resid. Dev F 41 198.232 1 15.488 40 182.744 7.2379 1 96.027 39 86.717 44.8752 1 9.332 38 77.384 4.3612 1 0.451 37 76.934 0.2107 Pr(>F)

NULL N_input 0.01065 Transectpos 7.164e-08 Noncrop_area 0.04370 N_input:Transectpos 0.64891

A kvázi modell

I

I I

I

A regressziós együtthatók ugyanazok és szignifikánsak maradtak. A standard hibák sokkal nagyobbak. Ha túlszóródásos modelleket hasonlítunk össze, akkor "Chi"-teszt helyett "F"-tesztet használunk. A túlszóródás egy lehetséges oka: az egyedek nem egymástól függetlenül, hanem csoportosan, klasztereződve jelennek meg.

Additív Modellek I

I

Hasznosak, ha a modellépítés folyamán nemlinearitást tapasztalunk. A lineáris modell nemparaméteres verziója. y = β0 +

p X

fj (Xj ) + ,

j=1

I I I

I I

ahol fj -k ”sima” függvények. Flexibilisebb modellek. Interpretálhatóak, hiszen fj -ket ki lehet rajzoltatni. Szimultán talál a magyarázó változóknak jó transzformációkat. Az R-ben legalább három csomag van rá (gam, mgcv, gss). Most a gam-ot használjuk.

Lineáris modell a gyomborítottságra

> mod5.lm <- lm(Weedcover ~ N_input + Transectpos + Noncrop_area > summary(mod5.lm) Call: lm(formula = Weedcover ~ N_input + Transectpos + Noncrop_area) Residuals: Min 1Q Median 3Q Max -25.67 -13.31 -6.19 11.12 42.08 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 30.733211 8.002910 3.840 0.000452 N_input -0.001777 0.036330 -0.049 0.961248 Transectpos1 -14.509641 5.369530 -2.702 0.010238 Noncrop_area 0.037191 0.159098 0.234 0.816425 Residual standard error: 17.4 on 38 degrees of freedom Multiple R-Squared: 0.1629, Adjusted R-squared: 0.09679 F-statistic: 2.464 on 3 and 38 DF, p-value: 0.0771

Reziduum–becsült érték ábra >plot(mod5.lm,1,pch=20) Residuals vs Fitted

40

39 ●

●

40 5●

● ● ●

20

● ● ●

Residuals

●

●

●

●

●

0

● ●

●

● ●

● ●

● ●

●

● ● ● ●

● ●

●

●

●

●

−20

●

●

● ●

●

●

●

20

25

Fitted values lm(Weedcover ~ N_input + Transectpos + Noncrop_area)

30

●

QQ-ábra >plot(mod5.lm,2,pch=20)

3

Normal Q−Q

39 ●

2

●

●

5

● ●

1

●

●

●

●

●

● ●

0

● ● ●

●

● ● ● ● ●

−1


● ●

40

● ●

●

●

●

●

●

●

●

●

●

● ●

●

● ● ●

●

●

−2

−1

0

1

Theoretical Quantiles lm(Weedcover ~ N_input + Transectpos + Noncrop_area)

2

Additív modell

> mod1.gam <- gam(Weedcover ~ lo(N_input) + lo(Noncrop_area) + + Transectpos)

Additív modell – summary

Call: gam(formula = Weedcover ~ lo(N_input) + lo(Noncrop_area) + Transec Deviance Residuals: Min 1Q Median 3Q Max -25.5007 -4.5381 -0.7272 4.6331 28.4130 (Dispersion Parameter for gaussian family taken to be 119.1652) Null Deviance: 13740.40 on 41 degrees of freedom Residual Deviance: 3559.746 on 29.8724 degrees of freedom AIC: 331.9167 Number of Local Scoring Iterations: 2 DF for Terms and F-values for Nonparametric Effects

(Intercept) lo(N_input) lo(Noncrop_area) Transectpos

Df Npar Df Npar F Pr(F) 1.0 1.0 5.0 10.1108 9.703e-06 1.0 3.1 4.5343 0.009068 1.0

Redukált modell a gyomborítottságra

> mod2.gam <- gam(Weedcover ~ lo(N_input) + lo(Noncrop_area))

Redukált modell a gyomborítottságra –summary Call: gam(formula = Weedcover ~ lo(N_input) + lo(Noncrop_area)) Deviance Residuals: Min 1Q Median 3Q Max -18.513 -8.883 -1.401 7.636 34.891 (Dispersion Parameter for gaussian family taken to be 178.6149) Null Deviance: 13740.40 on 41 degrees of freedom Residual Deviance: 5514.262 on 30.8724 degrees of freedom AIC: 348.2979 Number of Local Scoring Iterations: 2 DF for Terms and F-values for Nonparametric Effects Df Npar Df Npar F Pr(F) (Intercept) 1.0 lo(N_input) 1.0 5.0 6.6791 0.0002536 lo(Noncrop_area) 1.0 3.1 3.5089 0.0253447

A modellek összehasonlítása > anova(mod1.gam, mod2.gam, test = "F") Analysis of Deviance Table Model 1: Weedcover ~ lo(N_input) + lo(Noncrop_area) + Transectpos Model 2: Weedcover ~ lo(N_input) + lo(Noncrop_area) Resid. Df Resid. Dev Df Deviance F Pr(>F) 1 29.872 3559.7 2 30.872 5514.3 -1.000 -1954.5 16.402 0.0003349

R-négyzet= > 1 - 3560/13740 [1] 0.7409025

A lineáris modell R-négyzete: 0.16

A magyarázó változók transzformációi >par(mfrow=c(1,3),pty="s") plot(mod1.gam,residuals=T,se=T,pch=20) par(mfrow=c(1,1))

0

●

● ● ●

●

●

●

● ●

●

● ● ● ● ● ●

●

● ●●

● ●

●

● ●

● ●

●

● ●

●

●

●

● ●

20

●

● ●

● ●

10

●

● ● ●

●

●

●

●

●

−20

● ●

●

●

−20

●

●

●

● ● ●●

● ● ●

●

●

●●

0

50

100

150 N_input

200

250

−20

−30

●

●

●

0

10

20

30

40

Noncrop_area

●

● ●●

50

60

● ● ●

0

●

●

−10

●

● ● ●

●

−10

0

●

●

●


●

●

●

30

30 20

● ●

10

10

● ●

●

0

● ●

lo(Noncrop_area)

20

● ●

●

lo(N_input)

●

●

● ●

1 ●

●

−10

30

●

●

●● ●● ● ●● ●● ● ● ● ● ● ●● ●

●

Transectpos

Reziduum-becsült érték ábra

30

>plot(predict(mod1.gam),residuals(mod1.gam),pch=20)

20

●

●

10

●

●

●

●

● ●

● ●

0

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

● ●

−10

●

●

●

● ● ●

●

−20

residuals(mod1.gam)

●● ● ● ●

●

0

10

20

30 predict(mod1.gam)

40

50

Normalitás vizsgálat >qqnorm(residuals(mod1.gam),pch=20) 30

Normal Q−Q Plot

20

●

●

●

10

●

●

0

● ●

●

●

●

● ●

● ● ●

●

●

●

●

●

●

● ● ● ● ● ● ● ●

●

−10

●

●

●

● ● ●

●

−20

Sample Quantiles

● ● ● ●

●

−2

−1

0 Theoretical Quantiles

1

2

Általánosított Additív Modellek

I

I I

Az Additív Modellek általánosítása nem normális eloszlású függő változóra. Az additív modelléhez hasonló struktúra. A scale argumentum határozza meg, hogy a diszperziós paramétert becsülni kell-e (scale=-1).

Általánosított Additív modell az összfajszám függő változóra

> mod3.gam <- gam(Totspeciesnb ~ lo(N_input) + lo(Noncrop_area) + Transectpos, family = quasipoisson, scale = -1)

Summary Call: gam(formula = Totspeciesnb ~ lo(N_input) + lo(Noncrop_area) + Transectpos, family = quasipoisson, scale = -1) Deviance Residuals: Min 1Q Median 3Q Max -2.1274 -1.0111 -0.2104 0.8387 2.2575 (Dispersion Parameter for quasipoisson family taken to be 1.8426) Null Deviance: 198.2319 on 41 degrees of freedom Residual Deviance: 54.9228 on 29.8771 degrees of freedom AIC: NA Number of Local Scoring Iterations: 6 DF for Terms and F-values for Nonparametric Effects

(Intercept) lo(N_input) lo(Noncrop_area) Transectpos

Df Npar Df Npar F Pr(F) 1.0 1.0 5.0 2.1583 0.08562 1.0 3.1 1.8639 0.15525 1.0

A magyarázó változók transzformációi >par(mfrow=c(1,3),pty="s") plot(mod3.gam,residuals=T,se=T,pch=20) par(mfrow=c(1,1))

0

2

●

1

●

● ●

2

●

●

2

●

●

●

● ●

−1

●

●

● ●

●

●

●

● ●

● ●

N_input

200

250

1

●

● ●

●● ● ● ●

●

●

10

●●

●

●

● ●● ●

20

30

40

Noncrop_area

●

●● ●

● ● ●

●

0

● ●

●

●

150

●● ●

● ●

0

1

●

●

●

●

●

● ●

●

100

●

●

●

●

●

● ●

●

●

50

●

● ●

● ●

● ●

0

●

●

● ●

●

−1

●

●

● ●

● ● ●

● ● ● ● ●

−2

0

● ● ● ●

●

●

●

●

0

● ●

lo(Noncrop_area)

●

●

−1

●

−2

lo(N_input)

●

● ●

−2

1

● ● ●


●

●

● ● ●

50

60

Transectpos

●

Reziduum–becsült érték ábra >plot(predict(mod3.gam),residuals(mod3.gam),pch=20)

●

2

●

● ●

●

●

●

●

●

● ●

●

● ●

● ● ●

0

● ●

●

−1

●

●

●

●

●

●

●

●

●

●

● ● ●

● ● ● ●

−2

residuals(mod3.gam)

1

●

●

● ●

2.5

3.0 predict(mod3.gam)

3.5

●

Normalitás vizsgálat >qqnorm(residuals(mod3.gam),pch=20) Normal Q−Q Plot ●

2

●

● ●

●

1

● ● ●

●

●

●

●

● ● ● ●

0

● ● ● ● ● ● ●

●

●

−1

● ● ●

●

● ● ●

●

●

● ● ●

−2

Sample Quantiles

● ●

●

● ●

−2

−1

0 Theoretical Quantiles

1

2

Fix és random hatások

I

”Jobb”-e az A populació, mint a B populáció? I

I

Két populációt szeretnénk összehasonlítani: fix hatás

Nagy-e a variabilitás a populációk között? I I I

Specifikus összehasonlítás nem érdekes. A populációk egy minta egy (meta-)populációból. A populációk közötti variabilitást akarjuk figyelembe venni: random hatás

ANOVA egy fix hatással

yi = µ0 + βj + i I I I I

j csoportok, βj -k különbözőek, βj -k függetlenek. βj -ket megbecsülhetjük.

A megfigyelések közötti korrelációk

I

Tekintsük βj -ket valószínűségi változóknak. Legyen β ∼ N (0, σβ2 ).

I

Így kétféle varianciánk van: I I

σβ2 -csoportok közötti σ 2 -megfigyelések közötti

I

Egy csoporton belüli megfigyelések közötti korreláció:

I

σβ2 /(σβ2 + σ 2 ) Intraclass correlation

Miért használjunk random hatásokat?

I

A random modellek megmondhatják, hogy: I I

I

A különböző forrásokból mennyi variancia származik. Milyen megfigyelések valószínűek egy új populációban?

Segíthetik a becslést: I

Két hasonló populációban a megfigyelések is hasonlóak.

Az adatok véletlen variabilitása I I

A megfigyelési egységek (subject) gyakran csoportosulnak. Ez a csoportosulás származhat abból, hogy I

I

I

I

I

az adatoknak hierarhikus vagy beágyazott (nested) struktúrája van. Pl. a megfigyelési egységeket egy szervezeten belüli csoportokból választjuk, és feltételezzük, hogy a csoporthoz tartozás befolyásolhatja a megfigyeléseket. Ismételt méresek (repeated measures). Egy-egy megfigyelési egységen többször végzünk megfigyelést. Ha a megfigyelések időbeli ismétlések, akkor longitudinális adatokról beszélünk. Bizonyos esetekben a random variabilitás mindkét fajtája jelen van. Közös tulajdonság: korrelált megfigylések adott csoporton belül. A függetlenséget feltételező elemzések nem megfelelőek.

Kevert modellek (Mixed effects models) A kevert modellek fix és random hatásokat is tartalmaznak: fix hatások vagy az egész populációra, vagy annak bizonyos részeire vonatkoznak (kezelés vs. kontroll, hím vs. nőstény stb); random hatások bizonyos megfigyelési egységekre (subjects) vonatkoznak. I

I

I

A lineáris kevert modell (lme) egy olyan lineáris modell, amely mind fix, mind random hatásokat tartalmaz. Az általánosított lineáris kevert modell (glmm) egy olyan általánosított lineáris modell, amely mind fix, mind random hatásokat tartalmaz. A nemlineáris kevert modell (nlme) egy olyan nemlineáris modell, amely mind fix, mind random hatásokat tartalmaz.

Kevert modellek illesztése

I

I

Restricted Maximum Likelihood (REML) módszerrel történik az illesztés. A modell: Y = Xβ + Zγ + Fix hatások + Random hatások + Hibatag

Random hatások a földhasználati példában

I

I

I

I I

A különböző gazdálkodók (Farmer) és transzekt párok hatása nem lett figyelembe véve. Ha ezeket fixeknek tekintjük, akkor becsülhetjük a hatásukat a várható értékben. Ha random hatásoknak tekintjük, akkor feltételezzük, hogy a várható értékük 0 . A hatásuk a variancia struktúrában jelentkezik. Most csak a Farmer faktort vesszük figyelembe.

40 30 20 10 0

Weedcover

50

60

70

>plot(Weedcover∼Farmer,pch=20)

AG

ET

NL Farmer

PP

SD

Lineáris modell

> mod1.lm <- lm(Weedcover ~ N_input * Transectpos + + Noncrop_area + Farmer)

Lineáris modell – summary > summary(mod1.lm) Call: lm(formula = Weedcover ~ N_input * Transectpos + Noncrop_area + Farmer) Residuals: Min 1Q Median 3Q Max -27.109 -7.041 -0.884 5.508 40.086 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 32.530954 14.569454 2.233 0.03246 N_input -0.091177 0.075699 -1.204 0.23698 Transectpos1 -17.866460 6.021901 -2.967 0.00556 Noncrop_area 0.009736 0.156296 0.062 0.95070 FarmerET -6.453820 8.631873 -0.748 0.45995 FarmerNL 24.825768 13.035739 1.904 0.06560 FarmerPP -2.393432 8.919498 -0.268 0.79011 FarmerSD 22.186274 8.247480 2.690 0.01112 N_input:Transectpos1 0.034600 0.048059 0.720 0.47662 Residual standard error: 12.41 on 33 degrees of freedom Multiple R-Squared: 0.6303, Adjusted R-squared: 0.5407

Anova tábla

> anova(mod1.lm) Analysis of Variance Table Response: Weedcover Df Sum Sq Mean Sq F value Pr(>F) N_input 1 6.5 6.5 0.0425 0.8379581 Transectpos 1 2214.9 2214.9 14.3879 0.0006028 Noncrop_area 1 16.5 16.5 0.1074 0.7451385 Farmer 4 6342.6 1585.7 10.3004 1.587e-05 N_input:Transectpos 1 79.8 79.8 0.5183 0.4766213 Residuals 33 5080.0 153.9

Reziduum vs. becsült érték ábra >plot(mod1.lm,1,pch=20)

40

Residuals vs Fitted

39 ●

40

● ● ●

● ●

● ●

●

●

0

● ● ● ●

● ●

●

● ● ●

●

●

●

● ●● ●

●

● ●●

● ●

●

●

●

● ●

●

●

−20

Residuals

20

●

37 ●

10

20

30

40

Fitted values lm(Weedcover ~ N_input * Transectpos + Noncrop_area + Farmer)

50

Normalitás vizsgálat (QQ-plot) >plot(mod1.lm,2,pch=20)

4

Normal Q−Q

3

39 ●

2 1

●

●

●

●

●

●

●

●

●

●

0

● ●

● ● ● ● ● ●

● ● ● ● ●

−1

● ● ● ● ● ●

●

●

●

●

●

●

●

●

−2

●

●

●

37

−3


40 ●

−2

−1

0

1

Theoretical Quantiles lm(Weedcover ~ N_input * Transectpos + Noncrop_area + Farmer)

2

Szórás–becsült érték ábra >plot(mod1.lm,3,pch=20)

2.0

Scale−Location 39 ●

37 ● 40

● ●

●

1.0

● ● ●

● ● ● ●● ● ● ●

●

●

●

●

● ●

●

●

● ●

●

0.5

● ● ●

● ● ● ●

●

●

●

●

● ● ●

0.0


1.5

●

10

20

30

40

Fitted values lm(Weedcover ~ N_input * Transectpos + Noncrop_area + Farmer)

50

Problémák

I I

Enyhén nem állandó variancia. Nem független megfigyelések. (Transzekt párok, ugyanahhoz a gazdálkodóhoz tartozó területek.)

Kevert modellek az R-ben

Főbb csomagok: I

nlme:

I

Matrix, lme4:

I

Mindkettő szerzője: Douglas Bates et al. Különböző számítási módok, különböző random hatás formulák. Valószínűleg az lmer() hatékonyabb, de lényegesen nehezebb használni.

I

I

lme() lmer()

Kevert modell

> library(nlme) > mod1.lme <- lme(Weedcover ~ N_input * Transectpos + + Noncrop_area, random = ~1 | Farmer) > library(lme4) > mod1.lmer <- lmer(Weedcover ~ N_input * Transectpos + + Noncrop_area + (1 | Farmer)) random=∼1|Farmer : Egyetlen random hatás minden csoportnak: Farmer

Az lme modell > mod1.lme Linear mixed-effects model fit by REML Data: NULL Log-restricted-likelihood: -168.8030 Fixed: Weedcover ~ N_input * Transectpos + Noncrop_area (Intercept) N_input 34.68286414 -0.05117730 Transectpos1 Noncrop_area -17.86335610 0.04170213 N_input:Transectpos1 0.03469429 Random effects: Formula: ~1 | Farmer (Intercept) Residual StdDev: 13.39403 12.36867 Number of Observations: 42 Number of Groups: 5

Az lmer modell > mod1.lmer

Linear mixed-effects model fit by REML Formula: Weedcover ~ N_input * Transectpos + Noncrop_area + (1 | Farmer) AIC BIC logLik MLdeviance REMLdeviance 349.6 360 -168.8 336.8 337.6 Random effects: Groups Name Variance Std.Dev. Farmer (Intercept) 179.40 13.394 Residual 152.98 12.369 number of obs: 42, groups: Farmer, 5 Fixed effects: Estimate Std. Error t value (Intercept) 34.68293 11.50429 3.0148 N_input -0.05118 0.05951 -0.8600 Transectpos1 -17.86336 6.00315 -2.9757 Noncrop_area 0.04170 0.14433 0.2889 N_input:Transectpos1 0.03469 0.04791 0.7242 Correlation of Fixed Effects: (Intr) N_inpt Trnsc1 Nncrp_ N_input -0.736 Transectps1 -0.262 0.312

Az lme modell összegzése > summary(mod1.lme) Linear mixed-effects model fit by REML Data: NULL AIC BIC logLik 351.6061 362.8825 -168.8030 Random effects: Formula: ~1 | Farmer (Intercept) Residual StdDev: 13.39403 12.36867 Fixed effects: Weedcover ~ N_input * Transectpos + Noncrop_area Value Std.Error DF t-value (Intercept) 34.68286 11.504234 33 3.0147912 N_input -0.05118 0.059510 33 -0.8599827 Transectpos1 -17.86336 6.003157 33 -2.9756605 Noncrop_area 0.04170 0.144327 33 0.2889415 N_input:Transectpos1 0.03469 0.047909 33 0.7241658 p-value (Intercept) 0.0049 N_input 0.3960 Transectpos1 0.0054 Noncrop_area 0.7744

Az lmer modell összegzése > summary(mod1.lmer)

Linear mixed-effects model fit by REML Formula: Weedcover ~ N_input * Transectpos + Noncrop_area + (1 | Farmer) AIC BIC logLik MLdeviance REMLdeviance 349.6 360 -168.8 336.8 337.6 Random effects: Groups Name Variance Std.Dev. Farmer (Intercept) 179.40 13.394 Residual 152.98 12.369 number of obs: 42, groups: Farmer, 5 Fixed effects: Estimate Std. Error t value (Intercept) 34.68293 11.50429 3.0148 N_input -0.05118 0.05951 -0.8600 Transectpos1 -17.86336 6.00315 -2.9757 Noncrop_area 0.04170 0.14433 0.2889 N_input:Transectpos1 0.03469 0.04791 0.7242 Correlation of Fixed Effects: (Intr) N_inpt Trnsc1 Nncrp_ N_input -0.736 Transectps1 -0.262 0.312

Anova tábla A fix hatások feltételes anova teszttel tesztelhetők az lme modell esetén. > anova(mod1.lme) (Intercept) N_input Transectpos Noncrop_area N_input:Transectpos

numDF denDF F-value p-value 1 33 13.526606 0.0008 1 33 0.738573 0.3963 1 33 14.477856 0.0006 1 33 0.079834 0.7793 1 33 0.524416 0.4741

> anova(mod1.lmer) Analysis of Variance Table Df Sum Sq Mean Sq N_input 1 112.99 112.99 Transectpos 1 2214.88 2214.88 Noncrop_area 1 12.21 12.21 N_input:Transectpos 1 80.23 80.23

Redukált lme modell > mod2.lme <- lme(Weedcover ~ Transectpos, + random = ~1 | Farmer) > mod2.lme Linear mixed-effects model fit by REML Data: NULL Log-restricted-likelihood: -164.1596 Fixed: Weedcover ~ Transectpos (Intercept) Transectpos1 30.53005 -14.52381 Random effects: Formula: ~1 | Farmer (Intercept) Residual StdDev: 11.68188 12.23175 Number of Observations: 42 Number of Groups: 5

Redukált lmer modell > mod2.lmer <- lmer(Weedcover ~ Transectpos + + (1 | Farmer)) > mod2.lmer Linear mixed-effects model fit by REML Formula: Weedcover ~ Transectpos + (1 | Farmer) AIC BIC logLik MLdeviance REMLdeviance 334.3 339.5 -164.2 338.0 328.3 Random effects: Groups Name Variance Std.Dev. Farmer (Intercept) 136.47 11.682 Residual 149.62 12.232 number of obs: 42, groups: Farmer, 5 Fixed effects: Estimate Std. Error t value (Intercept) 30.530 5.899 5.176 Transectpos1 -14.524 3.775 -3.848 Correlation of Fixed Effects: (Intr) Transectps1 -0.320

Modellek összehasonlítása

> anova(mod1.lmer, mod2.lmer)

Data: Models: mod2.lmer: Weedcover ~ Transectpos + (1 | Farmer) mod1.lmer: Weedcover ~ N_input * Transectpos + Noncrop_area + (1 | Farme Df AIC BIC logLik Chisq Chi Df mod2.lmer 3 344.04 349.26 -169.02 mod1.lmer 6 348.82 359.24 -168.41 1.2253 3 Pr(>Chisq) mod2.lmer mod1.lmer 0.747

Az illesztett modell diagnosztikája

I I I

Jó kiindulás a modell fix verziójának ellenőrzése. Reziduumok ellenőrzése. Random faktor ellenőrzése.

>plot(mod2.lme) ●

3

●


2

●

●

1 ●

● ●

●

● ● ●

0

● ● ● ●

●

●

●

●

● ●

●

● ●

●

● ●

● ●

●

● ● ●

● ●

●

−1

● ●

−2 ●

10

20

30

Fitted values

40

A hibatag normalitásának ellenőrzése >qqnorm(mod2.lme) ●

2 ● ● ● ●

Quantiles of standard normal

● ●

1

● ● ●

0 ● ● ● ● ●

● ● ●

●

● ●

● ● ● ● ● ● ●

●

●

● ●

● ● ● ●

−1 ● ● ● ● ●

−2 ●

−2

−1

0


1

2

3

A random hatások becslése

> random.effects(mod2.lme) AG ET NL PP SD

(Intercept) -5.017905 -8.797690 4.564477 -7.836253 17.087371

A random hatás normalitásának ellenőrzése >qqnorm(mod2.lme, ranef(.)) (Intercept) ●


1.0

0.5

●

0.0

●

−0.5

●

−1.0 ●

−10

−5

0

5

Random effects

10

15

Problémák

I I I I

Van három kiugró érték. Nem állandó variancia. Nem normális reziduumok. Nem normális random hatások.

Variancia függvények I

I

I

A csoporton belüli hiba variancia struktúra valamilyen kovariáns függvényében. varFunc osztályok: különböző variancia függvények megadása: Argumentum: egyoldali formula, amely meghatározza a I I

I I

I

a kovariánst; opcionális rétegző változót: a rétegző változó minden szintjéhez más variancia paraméter tartozik. e.g. form = age| sex form = fitted(.) - a kovariáns = becsült értékek.

Standard varFunc osztályok: I I

I I

varFixed - lineáris; varIdent - különböző variancia a rétegző változó minden szintjéhez; varPower - hatvány variancia függvény. ...

> (mod3.lme <- lme(Weedcover ~ N_input + Transectpos, + weights = varIdent(form = ~1 | Farmer), + random = ~1 | Farmer)) Linear mixed-effects model fit by REML Data: NULL Log-restricted-likelihood: -152.6822 Fixed: Weedcover ~ N_input + Transectpos (Intercept) N_input Transectpos1 35.89923846 -0.05237175 -16.05589320 Random effects: Formula: ~1 | Farmer (Intercept) Residual StdDev: 12.92851 7.388345 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | Farmer Parameter estimates: AG ET NL PP SD 1.0000000 0.7059728 3.6631211 0.9700587 1.1778046 Number of Observations: 42 Number of Groups: 5

summary(mod3.lme) Linear mixed-effects model fit by REML Data: NULL AIC BIC logLik 323.3643 338.3364 -152.6822 Random effects: Formula: ~1 | Farmer (Intercept) Residual StdDev: 12.92851 7.388345 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | Farmer Parameter estimates: AG ET NL PP SD 1.0000000 0.7059728 3.6631211 0.9700587 1.1778046 Fixed effects: Weedcover ~ N_input + Transectpos Value Std.Error DF t-value p-value (Intercept) 35.89924 7.221607 35 4.971087 0.000 N_input -0.05237 0.036305 35 -1.442558 0.158 Transectpos1 -16.05589 2.189421 35 -7.333398 0.000 Correlation: (Intr) N_inpt N_input -0.503

Standardizált reziduumok vs. becsült értékek >plot(mod3.lme) ● ●

●

● ● ●

1

● ●

●

●

●


●

●

●

● ●

●

●

● ● ●

●

0 ●

●

●

● ● ●

●

●

●

●

● ● ●

●

−1 ●

● ● ●

−2 10

20

30

Fitted values

40

50

A hibatag normalitásának ellenőrzése >qqnorm(mod3.lme) ●

2 ● ● ● ●


● ● ●

1 ● ● ●

0

● ● ●

● ● ●

● ●

●

● ●

●

●

●

●

● ●

● ●

●

●

● ● ●

−1

● ● ● ● ● ●

−2 ●

−2

−1

0


1

A random hatás normalitásának ellenőrzése >qqnorm(mod3.lme,∼ranef(.)) (Intercept) ●


1.0

0.5

●

0.0

●

−0.5

●

−1.0 ●

−10

−5

0

5

Random effects

10

15

Következtetések

I I I

I

Csak a transzekt pozíciójának van szignifikáns hatása. A Farmer és N input hatása nem elválasztható. Több N input szintre lenne szükségünk minden Farmer esetén ahhoz, hogy szétválaszthassuk a hatásokat. Az additív modell lo(N input) és lo(Noncrop area) tényezője szignifikáns. A Farmer hatást nem vettük be a modellbe!

Általánosított lineáris kevert modell (Generalized Linear Mixed Model)

I

I I

Lineáris kevert modell nemnormális eloszlású függő változóval. Több függvény is van ilyen modellek illesztésére: glmmPQL (MASS package) I I

I I

Az lme és glm függvényekre épül. Ugyanolyan a struktúrája. A túlszóródás (overdispersion) is modellezhető.

glmmML (glmmML

package) lmer (Matrix package)

Általánosított lineáris kevert modell

> library(MASS) > mod1.glmm <- glmmPQL(Totspeciesnb ~ N_input + + Noncrop_area + Transectpos, random = ~1 | + Farmer, family = poisson, data = land) > mod1.glmm.lmer <- lmer(Totspeciesnb ~ N_input + + Noncrop_area + Transectpos + (1 | Farmer), + family = poisson, data = land)

> summary(mod1.glmm) Linear mixed-effects model fit by maximum likelihood Data: land AIC BIC logLik NA NA NA Random effects: Formula: ~1 | Farmer (Intercept) Residual StdDev: 5.462416e-06 1.375773 Variance function: Structure: fixed weights Formula: ~invwt Fixed effects: Totspeciesnb ~ N_input + Noncrop_area + Transectpos Value Std.Error DF t-value p-value (Intercept) 3.360393 0.13574646 34 24.754919 0.0000 N_input -0.001129 0.00065716 34 -1.718394 0.0948 Noncrop_area 0.005739 0.00271202 34 2.116053 0.0417 Transectpos1 -0.632091 0.09592256 34 -6.589594 0.0000 Correlation: (Intr) N_inpt Nncrp_ N_input -0.662 Noncrop_area -0.820 0.361 Transectpos1 -0.252 0.004 0.009

> summary(mod1.glmm.lmer)

Generalized linear mixed model fit using Laplace Formula: Totspeciesnb ~ N_input + Noncrop_area + Transectpos + (1 | Farm Data: land Family: poisson(log link) AIC BIC logLik deviance 87.14 95.83 -38.57 77.14 Random effects: Groups Name Variance Std.Dev. Farmer (Intercept) 0.002051 0.045288 number of obs: 42, groups: Farmer, 5 Estimated scale (compare to 1 ) 1.354733 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.3634085 0.1050857 32.01 < 2e-16 N_input -0.0010584 0.0005101 -2.07 0.03800 Noncrop_area 0.0053966 0.0019799 2.73 0.00642 Transectpos1 -0.6329749 0.0663255 -9.54 < 2e-16 Correlation of Fixed Effects: (Intr) N_inpt Nncrp_ N_input -0.697 Noncrop_are -0.804 0.389

Boxplot

30

40

50

>boxplot(Totspeciesnb∼Farmer)

10

20

●

AG

ET

NL

PP

SD

Redukált modell > (mod2.glmm <- glmmPQL(Totspeciesnb ~ Noncrop_area + + Transectpos, random = ~1 | Farmer, family = poisson)) Linear mixed-effects model fit by maximum likelihood Data: land Log-likelihood: NA Fixed: Totspeciesnb ~ Noncrop_area + Transectpos (Intercept) Noncrop_area Transectpos1 3.213777951 0.006774636 -0.631679342 Random effects: Formula: ~1 | Farmer (Intercept) Residual StdDev: 0.07566479 1.374297 Variance function: Structure: fixed weights Formula: ~invwt Number of Observations: 42 Number of Groups: 5

> summary(mod2.glmm) Linear mixed-effects model fit by maximum likelihood Data: land AIC BIC logLik NA NA NA Random effects: Formula: ~1 | Farmer (Intercept) Residual StdDev: 0.07566479 1.374297 Variance function: Structure: fixed weights Formula: ~invwt Fixed effects: Totspeciesnb ~ Noncrop_area + Transectpos Value Std.Error DF t-value p-value (Intercept) 3.213778 0.11003082 35 29.207979 0.0000 Noncrop_area 0.006775 0.00263462 35 2.571393 0.0145 Transectpos1 -0.631679 0.09454454 35 -6.681288 0.0000 Correlation: (Intr) Nncrp_ Noncrop_area -0.796 Transectpos1 -0.304 0.008 Standardized Within-Group Residuals:

>plot(mod2.glmm) ● ●

2

● ● ●


●

1

●

●

● ●

●

●

● ●

● ●

●

● ●

●

0 ● ● ●

●

● ●

●

● ●

●

● ●

−1

●

●

●

●

● ●

● ●

●

2.6

2.8

3.0

3.2

Fitted values

●

3.4

3.6

A hibatag normalitásának ellenőrzése >qqnorm(mod2.glmm) ●

2 ● ● ● ●


● ●

1

● ● ● ●

0

●

● ●

● ● ●

● ● ●

●

● ●

●

●

● ●

● ●

●

●

● ● ● ● ●

−1 ● ● ● ● ●

−2 ●

−1

0


1

2

A random hatás normalitásának ellenőrzése >qqnorm(mod2.glmm,∼ranef(.)) (Intercept) ●


1.0

0.5

●

0.0

●

−0.5

●

−1.0 ●

−0.05

0.00

Random effects

0.05

A hernyók méretbeli különbségei

I

I

I

Longitudinális adatokra példát szolgáltatnak egy olyan tömegnövekedést vizsgáló laboratóriumai kísérletből származó mérési eredmények, amelyben farkasalamlepke (Zerynthia polyxena) hernyóit különböző hőmérsékleten tartották. A teljes adathalmaz itt bemutatott részében 20 hernyó adatai szerepelnek (subjects). A hernyók testtömegét 4 naponta mérték a kikeléstől a bebábozódásig. 3 különböző hőmérsékleti kezelést alkalmaztak.

"Széles" formátumú táblázat

> mass[1:5, ] 1 2 3 4 5

BOX 3 4 7 8 9

TEMPR cooled cooled cooled cooled cooled

LM0 0.006 0.014 0.026 0.026 0.012

LM1 0.027 0.067 0.038 0.053 0.023

LM2 0.079 0.157 0.089 0.107 0.078

LM3 0.109 0.200 0.110 0.137 0.122

LM4 0.223 0.342 0.225 0.239 0.164

LM5 0.279 0.287 0.276 0.329 0.280

LM6 0.371 0.310 0.344 0.375 0.397

LM7 0.288 0.290 0.352 0.392 0.333

LM8 0.290 0.280 0.300 0.322 0.320

A széles típusú elrendezés átrendezése hosszúvá ("wide" → "long") > > > + + >

mass$BOX <- factor(mass$BOX) nobs <- nrow(mass) mass_long <- reshape(mass, idvar = "BOX", varying = list(names(mass[, 3:11])), direction = "long", , v.names = "LM") mass_long[1:10, ]

3.1 4.1 7.1 8.1 9.1 10.1 13.1 14.1 16.1 18.1

BOX 3 4 7 8 9 10 13 14 16 18

TEMPR time LM cooled 1 0.006 cooled 1 0.014 cooled 1 0.026 cooled 1 0.026 cooled 1 0.012 cooled 1 0.025 room 1 0.004 room 1 0.015 room 1 0.007 room 1 0.025

groupedData objektumok

I I

I

data.frame objektum kiterjesztése. Egy formulába foglalja egy függő változót, egy kovariánst és egy csoportosító tényezőt : függő ∼ kovariáns|csoport. Meg lehet külső faktorokat is adni: I

I I

a megfigyelési egységeken alkalmazott kezelés.

Cimkék is definiálhatók. ...

Egy groupedData objektum megadása

> library(nlme) > massGrouped <- groupedData(LM ~ time | BOX, outer = ~TEMPR, + data = mass_long, + labels = list(x = "Mérés", y = "Hernyótömeg (g)"))

>plot(massGrouped) 2

4

30 0.4

●

2

● ● ●

●

0.1 ●

●

●

●

23

●

●

●

25

● ● ●

●

●

● ● ●

0.0

● ●

●

2

4

6

8

●

6

●

●

●

●

● ● ●

2

Mérés

● ● ● ●

● ● ●

8

● ●

● ●

● ●

4

0.0

●

●

●

● ●

0.1

●

21

●

● ●

●

●

● ●

9 ●

●

2

●

● ●

●

●

0.2

●

●

8 ●

●

●

●

● ●

● ●

● ●

● ●

0.3

●

●

●

3

●

●

0.4 ● ●

●

●

●

7

●

0.1

●

●

●

● ● ●

● ●

● ● ●

●

● ●

0.4 ●

24 ●

●

10 ●

31 ● ● ●

●

● ●

4

●

●

●

●

●

●

28

● ●

● ● ●

● ●

●

●

●

●

● ●

●

●

●

● ● ●

● ● ●

●

22

●

● ●

8

●

● ●

●

●

26

6

18

●

●

●

4

● ●

●

●

● ●

0.0

0.2

2

13

●

●

●

0.3

8

●

● ●

●

● ●

6

20 ●

●

●

4

14

●

0.2

Hernyótömeg (g)

8

●

●

0.3

6

16 ●

4

6

8

● ● ●

● ●

2

4

6

8

Az ábra tulajdonságai

I I I

I

Ábra sorozat. Az egyedi ábrákat "panel"-eknek nevezzük. Minden panelnek egyedi tengelyei vannak. Csoportos adatok esetén a függő változónak egy elsődleges kovariánstól való csoportonkénti függésének tulajdonságai Trellis grafika (library(lattice)) segítségével könnyen áttekinthetők. Ugyanakkor a csoportok is összehasonlíthatók.

3 ábra a három hőmérsékleti kezelésre >plot(massGrouped,outer=T,layout=c(3,1),key=F)

2

4

cooled 0.4

●

●

Hernyótömeg (g)

0.3 ●

● ● ● ● ● ● ● ● ●

0.1

0.0

● ● ●

● ● ● ● ● ●

2

● ● ● ●

● ● ●

● ●

● ●

● ● ● ●

●

●

● ● ●

● ● ● ● ●

● ● ● ● ●

●

●

● ● ●

● ● ● ● ● ●

4

room ●

●

● ●

8

●

● ●

0.2

6

heated

6

● ● ● ● ●

●

●

● ●

●

● ● ●

● ●

● ● ● ● ● ● ● ●

● ● ●

● ● ● ● ●

● ● ●

● ●

● ●

●

● ● ● ●

●

● ● ●

●

● ●

● ● ● ● ●

● ● ● ●

● ● ● ● ●

●

8

●

● ●

2

Mérés

● ● ●

● ● ●

4

6

8

● ● ● ● ●

Az adatok tulajdonságai I

Szigmoid alak. A növekedés logisztikus függvénnyel modellezzük. y(x) =

I I I I

I

Φ1 1 + exp ΦΦ2 −x 3

az alsó aszimptota 0; Φ1 - felső aszimptota; Φ2 - inflekciós pont; Φ3 - skálázási paraméter (az inflexiós pont és a Φ1 /(1 + e−1 ) ≈ 0.73Φ1 távolsága az x tengelyen). A hernyómérettel növekvő variabilitás. Csoporton belüli variancia modell.

A logisztikus függvény

Lineáris modell

> mod1.lme <- lme(LM ~ TEMPR, weights = varPower(), + data = massGrouped) > anova(mod1.lme) (Intercept) TEMPR

numDF denDF F-value p-value 1 160 366.2286 <.0001 2 17 1.0027 0.3876

Reziduumok ábrája >plot(mod1.lme) ● ● ●

● ● ●

●

● ● ● ●

1

● ●


● ●●

● ●

● ● ● ● ● ● ● ● ● ● ● ● ● ●

● ● ● ● ● ●

●

●

● ●

● ● ● ●

● ●

● ●

0

● ●

● ● ● ● ●

●

●●

● ●

● ●

● ●

● ● ●

●● ●

● ●

● ●

●

●

● ●● ●

●

●

● ● ● ●

● ●

●●

● ● ●

● ●

−1 ● ●

● ●

●

● ● ●● ● ● ●

0.16

● ●

● ●

● ●

● ●

●

●

0.17

● ●

● ● ● ● ● ● ●

● ●

●

●

●

● ●

● ●

● ●

●

●

●

0.18

Fitted values

● ● ● ●

●

● ● ●

● ● ● ● ● ●

● ●

0.19

0.20

Normalitás vizsgálat qqnorm(mod1.lme) 3 ●

● ●


2

1

0

●● ● ● ● ● ● ●●● ●● ●● ●●● ● ● ●●● ●● ● ● ● ● ● ●

−1

● ● ●● ●

●●● ● ●●●● ● ● ● ●●●

●● ●●● ●

●

● ● ●●●

●

●● ●● ●● ●●

●●

● ● ●● ● ● ● ●● ●● ● ● ● ● ●● ● ●● ● ● ● ● ● ●● ●● ● ● ● ● ● ●● ● ● ● ●● ●● ●●● ● ● ●●●● ● ●●

● ● ● ●● ●● ● ● ● ● ● ● ●

● ● ●

−2

●

●

●

−3 −1

0


1

● ● ●● ●● ●●

● ● ●

● ● ●

selfStart

I

I

I

objektumok

A nemlineáris regressziós modellek egy osztálya kisegítő "inicializáló" függvénnyel. Tartalmazza a modell függvényt magát és a deriváltjainak kódját. Az iterációk megkezdése előtt az inicializáló függvényt hívja meg megfelelő kezdeti becslések készítéséhez.

Egyedi illesztések > mod1.nlsList <- nlsList(LM ~ SSlogis(time, Asym, xmid, + scal), data = massGrouped) > coef(mod1.nlsList) 4 10 7 3 8 9 21 23 26 22 28 25 31 24 30 16 14 20

Asym 0.3013417 0.2895020 0.3423090 0.3185156 0.3766248 0.3593805 0.2246673 0.2595412 0.3133955 0.2629832 0.2845404 0.3226826 0.3384038 0.3649638 0.3723890 0.2923923 0.2882872 0.2910224

xmid 3.019328 3.591827 4.405085 4.295591 4.276009 4.739875 5.616000 5.136916 4.285689 5.294313 3.691796 4.554836 4.271248 4.906504 3.774257 4.830018 4.243042 3.422182

scal 0.7395170 0.7846719 1.0376891 0.8268333 1.0694629 0.9874677 1.0037107 1.0256070 1.1490976 0.7014190 0.6496305 0.8339809 1.0323817 0.8919346 0.8482556 1.1779122 1.2234200 1.0152101

>plot(intervals(mod1.nlsList),col=as.numeric(TEMPR)) Asym 18

xmid

|

13

|

|

20

|

14

25

BOX

|

|

|

22

23

| |

|

|

|

3

|

7

|

|

4

|

0.25

|

|

|

| |

|

|

| |

0.30

|

|

|

0.35

|

0.40

|

3

|

|

|

| |

|

|

| |

| |

|

|

|

|

|

| |

|

7

| |

|

|

|

6

|

|

|

5

|

| |

|

|

|

|

|

4

| |

|

|

|

| |

|

|

|

|

|

|

|

| |

|

| |

8

0.20

|

|

| |

| |

|

|

0.15

|

| |

|

|

|

|

|

|

|

|

|

|

|

|

|

|

|

|

|

|

|

9

10

|

| |

| |

|

|

|

|

|

| |

|

|

|

|

|

|

|

|

|

|

|

|

|

26

|

|

|

|

|

| |

|

|

|

|

|

|

28

|

|

|

|

|

|

|

|

24

scal

|

| |

30

21

| |

|

31

| |

|

|

|

|

|

|

|

16

|

|

|

|

0.5

|

1.0

1.5

2.0

nlme modell > (mod1.nlme <- nlme(mod1.nlsList, na.action = na.omit, + control = list(tolerance = 0.001))) Nonlinear mixed-effects model fit by maximum likelihood Model: LM ~ SSlogis(time, Asym, xmid, scal) Data: massGrouped Log-likelihood: 335.0960 Fixed: list(Asym ~ 1, xmid ~ 1, scal ~ 1) Asym xmid scal 0.3106634 4.3041546 0.9230615 Random effects: Formula: list(Asym ~ 1, xmid ~ 1, scal ~ 1) Level: BOX Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr Asym 0.03472737 Asym xmid xmid 0.58421223 -0.344 scal 0.07711092 -0.029 0.948 Residual 0.03105154 Number of Observations: 180 Number of Groups: 20

Populaciós és csoporton belüli becslések >plot(augPred(mod1.nlme,level=0:1)) fixed 2

4

16

6

BOX

8

2

14

4

20

6

18

● ●

●

●

●

●

●

●

●

●

●

●

●

●

25

Hernyótömeg (g)

●

●

●

●

●

● ●

●

●

●

●

9

● ●

●

●

●

●

●

●

●

●

● ●

●

● ●

0.0

30

●

●

0.1 0.0

● ●

●

●

●

●

● ●

21

●

●

●

23

●

●

26

22 0.4

● ●

●

●

●

● ●

●

● ● ●

●

●

●

●

●

0.0

●

●

● ● ●

● ●

●

2

4

● ●

6

8

●

● ●

● ●

●

●

●

● ●

4

0.0 ●

●

6

Mérés

8

●

● ●

● ●

●

● ●

●

●

2

0.2 0.1

8

●

●

2

●

●

●

●

●

●

●

●

●

0.3 ●

●

●

●

●

●

●

3 ●

●

●

●

●

●

●

●

●

0.1

●

●

●

●

●

7

0.4

0.2

●

●

10

●

●

●

●

4 0.3

●

●

● ●

0.3

0.1

●

24

● ●

●

●

●

31

● ●

●

0.2 ●

●

●

0.2

0.4 ●

●

●

0.4 0.3

●

●

●

●

28

●

●

●

●

● ●

●

●

●

●

●

●

●

●

● ●

8

13

4

6

8

Standardizált reziduumok vs. becsült értékek >plot(mod1.nlme,pch=20,col=1) ●

●

● ●

●

2

● ●

●

● ● ●

●

●

●


●

●

●

●

●

●

●

1

● ●

● ● ●● ● ●● ● ● ● ●● ●●

●

●

●●

● ●● ●

●

●

●

● ●

●

● ● ● ●● ● ●● ●●● ● ● ● ● ●

●

● ●

● ●

● ●

● ●

●

● ●

● ●

●

● ●

● ● ●● ● ● ●

● ● ●

●

● ● ●

●●

●

●

● ●

●

●

●

●

● ●● ● ● ● ●

● ●

●

●

●

●

● ● ● ●

● ● ● ●

−2 ●

0.0

0.1

0.2

Fitted values

●

●

● ●

●

●

●● ●

● ●

●

●

●

●

●

●

●

−1

●

●

● ●

●

● ●

●

●

●

● ●

● ●

●

●

0

●

● ●

● ●

●

●

●

● ●

●

● ●

0.3

●

●

●

A reziduumok normalitása >qqnorm(mod1.nlme,pch=20,col=1) 3 ●

● ● ●


2

1

0

−1

● ● ● ● ●

−2

● ●

●

● ●● ● ●● ● ●● ● ●● ● ● ●● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ●● ● ● ●● ● ●● ●

● ●● ● ● ●● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●

●● ●●● ●● ● ●● ●● ● ●● ● ● ●● ● ● ●● ● ● ●● ● ● ●● ●●● ● ●●

● ●● ● ●● ● ●● ● ● ●

●● ● ●● ●● ●

●

●

●

●

● ● ●

●

●

−3 −2

−1

0


1

2

A random hatások normalitása >qqnorm(mod1.nlme,∼ranef(.),pch=20,col=1) Asym

xmid

2

scal

●

●

●

●

●


1

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

−2

●

●

●

−1

●

●

●

0

●

●

●

●

●

−0.05

●

0.00

0.05

−1.0

●

−0.5

0.0

0.5

Random effects

1.0

−0.15

−0.10

−0.05

0.00

0.05

0.10

> (mod2.nlme <- update(mod1.nlme, + weights = varPower(form = ~time))) Nonlinear mixed-effects model fit by maximum likelihood Model: LM ~ SSlogis(time, Asym, xmid, scal) Data: massGrouped Log-likelihood: 378.8348 Fixed: list(Asym ~ 1, xmid ~ 1, scal ~ 1) Asym xmid scal 0.3274611 4.4281031 1.1078004 Random effects: Formula: list(Asym ~ 1, xmid ~ 1, scal ~ 1) Level: BOX Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr Asym 0.028874955 Asym xmid xmid 0.640498532 -0.379 scal 0.129289287 0.261 0.439 Residual 0.004261739 Variance function: Structure: Power of variance covariate Formula: ~time Parameter estimates: power

Reziduumok ábrája >plot(mod2.nlme,pch=20,col=1) ●

2 ●

● ● ●● ●

●

● ●

● ●

●

●

●

● ●

● ●

●

1

● ●


● ● ● ●

●

●

● ●

●

●

●

● ●

●

●

●

●

0

● ● ●

●

●

● ●

● ●

● ● ●

●

●

●

●

●

●● ●●

●

●

●

● ●

●

●

●

● ●

●

●

●

●●

●

●

●

●● ● ● ● ● ●

● ● ● ●

●

●

● ●

●

−2 ●

●

0.0

0.1

0.2

Fitted values

●

●

●

● ●

●

●

● ● ●

●

● ●

●

●

●●

●

●

● ●

●

● ●

● ● ● ●

●

−1

●

●

●

● ●

● ●

●

●

● ●

● ●

●

●

● ●● ●● ● ● ● ●● ●

●

●●

●

●

●

●

● ● ●●

● ●

● ●

●

●

●● ●

●

●

●

● ● ● ●

●

0.3

●

A random hatások normalitása >qqnorm(mod2.nlme,∼ranef(.),pch=20,col=1) Asym

xmid

2

scal

●

●

●

●

●


1

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

−2

●

●

●

−1

●

●

●

0

●

●

●

●

●

−0.04

●

−0.02

0.00

0.02

0.04

●

−1.0

−0.5

0.0

0.5

Random effects

1.0

1.5

−0.2

−0.1

0.0

0.1

Modellek összehasonlítása

> anova(mod1.nlme, mod2.nlme) mod1.nlme mod2.nlme mod1.nlme mod2.nlme

Model df AIC BIC logLik Test L.Ratio 1 10 -650.1921 -618.2625 335.0960 2 11 -735.6695 -700.5470 378.8348 1 vs 2 87.47744 p-value <.0001

Fix hatás megadása a modellben

> modfix <- fixed.effects(mod2.nlme) > mod3.nlme <- update(mod2.nlme, fixed = Asym + xmid + + scal ~ TEMPR, start = c(modfix[1], 0, 0, modfix[2], + 0, 0, modfix[3], 0, 0)) > anova(mod3.nlme) Asym.(Intercept) Asym.TEMPR xmid.(Intercept) xmid.TEMPR scal.(Intercept) scal.TEMPR

numDF denDF F-value p-value 1 152 689.1728 <.0001 2 152 3.2536 0.0413 1 152 270.7553 <.0001 2 152 1.6871 0.1885 1 152 626.2964 <.0001 2 152 0.2687 0.7648

>plot(mod3.nlme,pch=20,col=1) ●

2

● ●

●

● ●

●

●

●

● ●

●

● ● ●

●

●●


●

●

●

●

●

●

0

●

●●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

● ●

●

●●

●

●

●

●

● ●

● ●

● ● ●

−2

●

●

●●

●

●

●

●

0.1

●

●

●

●

●

● ● ● ● ● ● ● ●●

●

● ●

●

●

● ●

●

●

●

● ●

● ● ●

● ●

●

0.0

● ● ●

●

●

●

●

●

● ●

●

● ●

●

●

−1

●● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●● ● ● ●

● ●

●

●

●

●

●

● ●

●

●

●

●

● ●

● ●

●

●

●

●

1

●

● ●

●

●

0.2

Fitted values

0.3

● ●

●

●

>qqnorm(mod3.nlme,pch=20,col=1) 3 ●

● ● ●


2

1

0

−1

● ● ●

−2

●

● ●

● ●● ●● ●● ● ● ●

● ●● ●● ● ● ● ●

●● ●

●●● ● ● ●● ● ●● ● ● ● ● ●●● ● ● ●● ● ● ●● ●● ● ● ● ● ● ●● ●●● ●● ● ●●● ●● ●●● ● ● ●●● ● ● ●●● ● ● ● ● ●● ●

●● ● ● ●● ●● ●● ● ● ●● ● ● ● ●● ● ● ●● ● ● ● ●●● ● ● ● ● ● ●● ● ● ●● ●●● ● ● ●●

●● ● ●●

● ●●

●

● ●

● ● ●

●

● ●

●

−3 −2

−1

0


1

2

Nemlineáris modellépítés

I I

I

I I

I

Iteratív modellépítési folyamat. Csoporton belüli modellépítés (nlsList) a megfigyelések minden csoportjára. Ezekből az illesztésekből meghatározzuk, hogy a paraméterek hogyan változnak a csoportok között, illetve hogy függnek egyéb kovariánsoktól stb. Bonyolultabb modellek illesztése a fentiek alapján. Ha a csoporton belüli modell nemlineáris, kezdő becslések szükségesek. Ismételjük meg az illesztést a fix tényező beillesztésével.

A lineáris modellektől a nemlineáris kevert modellekig R-ben

Recommend Documents