12. nodaļa. Vispārinātie lineārie modeļi
12.1 Teorētiskais pamatojums
Vēl jāuzraksta!
Plašāku materiālu par vispārinātajiem lineārajiem modeļiem var skatīt Dobson and Barnett (2008) un McCulloch, Searle, and Neuhaus (2008).
12.2 Dati
Divu GLM analīzes veidu (binārā loģistiskā regresija un puasona regresija) piemēram izmantoti četri datu faili. Pirmais fails nezales.txt satur informāciju par pētījumu, kurā analizēts nezāļu skaits parauglaukumos atkarībā no pielietotā augu aizsardzības līdzekļa - kontrole, līdzeklis 1 un līdzeklis 2 (kolonna grupa), kā arī augsnes pH šajos parauglaukumos.
nezales<-read.table(file="nezales.txt",header=T)
str(nezales)## 'data.frame': 300 obs. of 3 variables:
## $ grupa : Factor w/ 3 levels "kontrole","lidzeklis1",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ pH : num 6.26 6.8 6.79 6.52 6.08 6.21 6.47 6.74 6.71 6.11 ...
## $ skaits: int 2 10 5 6 5 4 2 2 5 3 ...head(nezales)## grupa pH skaits
## 1 kontrole 6.26 2
## 2 kontrole 6.80 10
## 3 kontrole 6.79 5
## 4 kontrole 6.52 6
## 5 kontrole 6.08 5
## 6 kontrole 6.21 4summary(nezales)## grupa pH skaits
## kontrole :100 Min. :5.310 Min. : 0.00
## lidzeklis1:100 1st Qu.:6.388 1st Qu.: 0.00
## lidzeklis2:100 Median :6.700 Median : 1.00
## Mean :6.699 Mean : 2.14
## 3rd Qu.:7.040 3rd Qu.: 3.25
## Max. :7.920 Max. :10.00Failos veids1.txt, veids2.txt un veids3.txt ir viena un tā paša pētījuma dati, kas parādīti trīs dažādos veidos. Šajā pētījumā ir novērtēts vai augs ir bijis bojāts vai nē (slimības ietekmēts) atkarībā no tā šķirnes un audzēšānas veida (tunelis vai lauks). Failā veids1 ir izejas dati, kur katram augam ir sava rindiņa un kolonnā bojajums ir atzīmē 0 vai 1 atkarībā no tā, vai augs nav bojāts, vai ir bojāts. Failā veids2 dati jau ir apkopoti tādā veidā, ka katrai sķirnes un audzēšanas veida kombinācijai ir norādīts bojāto un nebojāto augu skaits. Failā veids3 dati arī jau ir apkopoti, tikai šajā gadījumā ir kolonna ar procentuālo bojāto augu daudzumu (izteikts decimāldaļās) un kopējo augu skaitu katrā faktoru kombinācijā
veids1<-read.table(file="veids1.txt", header=TRUE,sep="\t",dec=".")
str(veids1)## 'data.frame': 175 obs. of 3 variables:
## $ skirne : Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
## $ veids : Factor w/ 2 levels "lauks","tunelis": 1 1 1 1 1 1 1 1 1 1 ...
## $ bojajums: int 1 1 1 1 0 1 0 1 0 1 ...head(veids1)## skirne veids bojajums
## 1 A lauks 1
## 2 A lauks 1
## 3 A lauks 1
## 4 A lauks 1
## 5 A lauks 0
## 6 A lauks 1veids2<-read.table(file="veids2.txt", header=TRUE,sep="\t",dec=".")
str(veids2)## 'data.frame': 6 obs. of 4 variables:
## $ skirne : Factor w/ 3 levels "A","B","C": 1 1 2 2 3 3
## $ veids : Factor w/ 2 levels "lauks","tunelis": 1 2 1 2 1 2
## $ bojats : int 17 10 14 16 21 2
## $ nebojats: int 3 17 11 24 13 27head(veids2)## skirne veids bojats nebojats
## 1 A lauks 17 3
## 2 A tunelis 10 17
## 3 B lauks 14 11
## 4 B tunelis 16 24
## 5 C lauks 21 13
## 6 C tunelis 2 27veids3<-read.table(file="veids3.txt",header=TRUE,sep="\t",dec=".")
str(veids3)## 'data.frame': 6 obs. of 4 variables:
## $ skirne : Factor w/ 3 levels "A","B","C": 1 1 2 2 3 3
## $ veids : Factor w/ 2 levels "lauks","tunelis": 1 2 1 2 1 2
## $ procents: num 0.85 0.37 0.56 0.4 0.618 ...
## $ skaits : int 20 27 25 40 34 29head(veids3)## skirne veids procents skaits
## 1 A lauks 0.85000000 20
## 2 A tunelis 0.37037037 27
## 3 B lauks 0.56000000 25
## 4 B tunelis 0.40000000 40
## 5 C lauks 0.61764706 34
## 6 C tunelis 0.06896552 2912.3 Puasona regresija
Gadījumos, kad regresors (atkarīgais mainīgais) ir skaita dati, turklāt to vērtības ir mazas un sadalījums ir izteikti asimetrisks (daudz mazo vērtību, veidojas Puasona sadalījums), izmantot lineāro regresiju nebūtu pareizi, jo tiktu pārkāpti vairāki nosacījumi, kas izvirzīti šai analīzei. Analizējot skaita datus kā alternatīvu var izmantot vispārinātos lineāros modeļus (GLM) norādot, ka datiem ir Puasona atlikuma struktūra. Šo GLM veidu sauc arī par Puasona regresiju. GLM analīzi programmā R veic ar funckiju glm(). Funkcijā jānorāda trīs komponentes: formula, kas jāpārbauda (regresents atkarībā no viena vai vairākiem regresoriem un to kombinācijas), atlikumu struktūras veids, kas šajā gadījumā ir family=poisson("log"), kā arī datu objekts, kurā atrodas analizējamie mainīgie. Arguments "log" iekavās pie poisson() parāda kāda saistības funkcija ir izmantota, lai pārietu no oriģinālās saistības starp mainīgajiem uz lineāro saistību. Analīzes rezultātus var apskatīt ar funckijām summary() un anova(). GLM modeļos var izmantot gan nepātraukti variējošus regresorus, gan arī kategorijas regresorus.
Piemērā apskatīts kā nezāļu skaitu ietekmē piederība pie kādas no pētījumu grupas un augsnes pH, kā arī šo divu faktoru kombinācija.
mod<-glm(skaits~grupa*pH,family=poisson("log"), data=nezales)
summary(mod) ##
## Call:
## glm(formula = skaits ~ grupa * pH, family = poisson("log"), data = nezales)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1527 -0.8192 -0.2107 0.6433 3.2340
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.51970 0.62685 2.424 0.0153 *
## grupalidzeklis1 -2.21275 1.48261 -1.492 0.1356
## grupalidzeklis2 -3.70175 2.41941 -1.530 0.1260
## pH 0.01156 0.09301 0.124 0.9011
## grupalidzeklis1:pH 0.11340 0.21933 0.517 0.6051
## grupalidzeklis2:pH 0.14455 0.36004 0.401 0.6881
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 873.24 on 299 degrees of freedom
## Residual deviance: 309.79 on 294 degrees of freedom
## AIC: 879.63
##
## Number of Fisher Scoring iterations: 5anova(mod,test="Chisq")## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: skaits
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 299 873.24
## grupa 2 562.83 297 310.41 <2e-16 ***
## pH 1 0.23 296 310.18 0.6350
## grupa:pH 2 0.39 294 309.79 0.8237
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Secinājumi: pēc analīzes kopsavilkuma tabulas var spriest, ka statistiski būtiska ir tikai Intercept vērtība (nezāļu skaits būtiski atšķiras no 0), tomēr pārējo koeficientu vērtības nav būtiskas. GLM analīzē, līdzīgi kā tas ir ANCOVA, ja ir iekļauts kvalitatīvs/kategorijas mainīgais (grupa), tad tas tiek pārkodēts un salīdzināšana tiek veikta references līmeni (grupu), kas šajā gadījumā ir kontrole. Koeficienti, kas parādās kopsavilkuma tabulā, ir koeficienti lineārājai vienādojumam. Lai atgrieztos pie sākotnējām skaita vērtībām, šis lineārais vienādojums ir jāizmanto kā pākāpe eksponentam - pretējā darbība naturālajam logaritmam. Ar funkciju `anova() var apskatīt kā katrs faktors un kombinācija kopumā ietekmē analizējamo pazīmi. Pēc šīs tabulas varam secināt, ka būtiska ietekme ir grupai, bet nav pH un grupas:pH kombinācijai.
Tā kā pH ietekme nebija būtiska, šajā gadījumā varam izveidot otru modeli, kurā kā faktors ir tikai piederība pie grupas.
mod2<-glm(skaits~grupa,family=poisson("log"),data=nezales)
summary(mod2)##
## Call:
## glm(formula = skaits ~ grupa, family = poisson("log"), data = nezales)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1432 -0.8000 -0.1522 0.7063 3.1686
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.59737 0.04499 35.50 <2e-16 ***
## grupalidzeklis1 -1.44895 0.10317 -14.04 <2e-16 ***
## grupalidzeklis2 -2.73680 0.18240 -15.01 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 873.24 on 299 degrees of freedom
## Residual deviance: 310.41 on 297 degrees of freedom
## AIC: 874.24
##
## Number of Fisher Scoring iterations: 5Secinājumi: Tā kā visi koeficienti ir būtiski, var secināt, ka gan grupā līdzeklis 1, gan līdzeklis 2 nezāļu skaits ir mazākas nekā kontroles grupā. Jāatceras, ka šie koeficienti nav interpretējami tiešā veidā kā atšķirības starp vidējiem nezāļu skaitiem, jo ir izmantota log saistības funckija. Sagaidāmais nezāļu skaits kontroles grupā attiecīgi ir \(e^1.597\), līdzeklis 1 grupā \(e^{1.597-1.449}\) un līdzeklis 2 grupā \(e^{1.597-2.737}\).
Tā kā tika ieveidots modelis, kurā bija divi mainīgie un to kombinācija, un modelis ar vienu mainīgo, tad šos modeļus var salīdzināt savā starpā, lai secinātu, vai modeļa vienkāršošana ir bijusi pamatota. Divu pakārtotu/saistītu modeļu salīdzināšanu var veikt ar funckijā AIC(), vai anova(), kurai kā arguments norādīts Chisq tests.
anova(mod,mod2,test="Chisq")## Analysis of Deviance Table
##
## Model 1: skaits ~ grupa * pH
## Model 2: skaits ~ grupa
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 294 309.80
## 2 297 310.41 -3 -0.61319 0.8934AIC(mod,mod2)## df AIC
## mod 6 879.6302
## mod2 3 874.2434Secinājumi: izmantojot hī kvadrāta testu, konstatēts, ka atšķirība starp modeļiem nav būtiska (p vērtība ir 0,89), attiecīgi mēs varam izvēlēties vienkāršāko modeli, kurā ir tikai viens regresors. Arī salīdzinot pēc AIC vērtībām, otraijam modelism ir mazāka AIC vērtība un tas uzskatāms par piemērotāku.
12.4 Binārā loģistiskā regresija
Binārās loģistiskās regresijas veikšanai izmanto to pašu funkciju glm(), tikai šajā gadījumā kā family= jānorāda binomial un saistības funkcija ir logit (ko var arī mainīt pret probit). Formāts, kādā veidā norāda regresentu (atkarīgo mainīgo) ir atkarīgs no tā, kādā veidā ir pieejami dati. Binārajā loģistiskajā regresijā pārbauda vai iespējamības (ka iestāsies konkrētais notikums (būs bojājums)) būtiski atšķiras, atkarībā no ietekmējošiem faktoriem (gan kvantitatīviem, gan kvalitatīviem).
Ja ir izejas dati, kuros atkarīgais mainīgais sastāv no 0 un 1, tad funkcijā glm() šī kolonna ir jānorāda, kā atkarīgais mainīgais (šadā veidā dati ir pieejami failā veids1).
mod1<-glm(bojajums~skirne+veids, family=binomial("logit"),data=veids1)
summary(mod1)##
## Call:
## glm(formula = bojajums ~ skirne + veids, family = binomial("logit"),
## data = veids1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7901 -0.9376 -0.6037 1.1195 1.8932
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3773 0.3978 3.462 0.000536 ***
## skirneB -0.4604 0.4203 -1.095 0.273354
## skirneC -1.2396 0.4446 -2.788 0.005304 **
## veidstunelis -1.7476 0.3521 -4.963 6.95e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 241.31 on 174 degrees of freedom
## Residual deviance: 208.73 on 171 degrees of freedom
## AIC: 216.73
##
## Number of Fisher Scoring iterations: 4Ja dati ir jau apkopoti un tajos ir kolonnas, kas parāda pozitīvo un negatīvo notikumu skaitu, kā failā veids2 kolonnas bojats un nebojats, tad kā atkarīgais mainīgais formulā ir jānorāda abas kolonnas, kas apvienotas ar funkciju cbind().
mod2<-glm(cbind(bojats,nebojats)~skirne+veids,family=binomial("logit"),data=veids2)
summary(mod2)##
## Call:
## glm(formula = cbind(bojats, nebojats) ~ skirne + veids, family = binomial("logit"),
## data = veids2)
##
## Deviance Residuals:
## 1 2 3 4 5 6
## 0.5944 -0.4049 -1.6402 1.2947 0.9796 -1.5666
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3773 0.3978 3.462 0.000536 ***
## skirneB -0.4604 0.4203 -1.095 0.273354
## skirneC -1.2396 0.4446 -2.788 0.005304 **
## veidstunelis -1.7476 0.3521 -4.963 6.95e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 40.8849 on 5 degrees of freedom
## Residual deviance: 8.2977 on 2 degrees of freedom
## AIC: 37.095
##
## Number of Fisher Scoring iterations: 4Trešā iespēja ir, ka dati apkopoti un tajos ir kolonna ar procentiem (decimāldaļās), kas parāda notiku iestāšanās biežumu (šajā gadījumā procentuālais bojāto augu daudzums). Lai analizētu šādus datus, kolonna ar procentiem jānorāda kā atkarīgais mainīgais, plus papildus jānorāda arguments weights= un kolonna, kurā ir oriģinālie kopējie skaita dati. To var nenorādīt gadījumos, ja skaita dati, no kuriem aprēķināti procenti ir vienādi starp visām faktoru kombinācijām (visās rindiņās).
mod3<-glm(procents~skirne+veids,weights=skaits,family=binomial("logit"),data=veids3)
summary(mod3)##
## Call:
## glm(formula = procents ~ skirne + veids, family = binomial("logit"),
## data = veids3, weights = skaits)
##
## Deviance Residuals:
## 1 2 3 4 5 6
## 0.5944 -0.4049 -1.6402 1.2947 0.9796 -1.5666
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3773 0.3978 3.462 0.000536 ***
## skirneB -0.4604 0.4203 -1.095 0.273354
## skirneC -1.2396 0.4446 -2.788 0.005304 **
## veidstunelis -1.7476 0.3521 -4.963 6.95e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 40.8849 on 5 degrees of freedom
## Residual deviance: 8.2977 on 2 degrees of freedom
## AIC: 37.095
##
## Number of Fisher Scoring iterations: 4Secinājumi: neatkarīgi no datu sagatavošanas veida, visi iegūtie koeficienti un tie būtiskumi ir vienādi. Pastāv būtiska atšķirība starp skirnes A (izmantota kā bāzes līmenis) un skirnes Cbojājumu iespējamību (šķirnei C tā ir mazāka nekā šķirnei A, ko parāda negatīvs koeficients pie šķirnes C), kā arī it būtiska atšķirība starp audzēšanas veidu - lauks un tunelis, turklāt tuneļa audzēšanas veidam bojājumu daudzums ir mazāks nekā laukam veidam. Jāatceras, ka arī šajā gadījumā šie ir lineārā vienādojuma koeficienti, kur y vērtības iegūtas izmantojot logit saistības funkciju no sākotnējiem bojājumu datiem.
Lai aprēķinātu prognozēto bojājumu iespējamību, var izmantot funkciju predict(), kurai kā arguments jānorāda modeļa objekts. Ja nenorāda papildus argumentus, tad aprēķinātās vērtības būs lineārajam vienādojumam. Prognozēto iespējamības vērtību iegūšanai jānorāda papildus arguments type="respone".
predict(mod3)## 1 2 3 4 5 6
## 1.3773246 -0.3702801 0.9168930 -0.8307117 0.1376954 -1.6099093predict(mod3,type="response")## 1 2 3 4 5 6
## 0.7985610 0.4084733 0.7144086 0.3034946 0.5343696 0.1666012Secinājumi: Izmantojot trešo modeli, iegūtas prognozētās bojājumu iespējamības vērtības katrai no sešām pārbaudītajām abu faktoru kombinācijām, piemēram, ja audzēšanas veids ir lauks un šķirne ir A (pirmā rindiņa failā veids3), tad prognozētais bojājumu iespējamība ir 0,798 jeb 79,8%.