Kordamine

1. Tehke PISA andmestiku alusel regressioonanalüüs, kus sõltuvaks tunnuseks on matemaatika testi tulemus (PVMATH) ja prediktoriteks enese-tõhusus teaduses (SCIEEFF; science self-efficacy) ja mina-pilt loodusteadustes (SCSCIE; science self-concept).

2. Arvutage ka standardiseeritud regressioonikordajad ja mudeli parameetrite usalduspiirid.

3. Kas mudeli jäägid jaotuvad normaalajaotusel vastavalt?

Logistiline regressioon

Kui lineaarse regressiooni puhul oli sõltuv tunnus numbriline, siis logistilise regressiooni puhul on sõltuv tunnus kategoriaalne. Logistiline mudel ennustab prediktorite väärtuste abil võimalust kuuluda mingisse kategooriasse. Binaarne logistiline regressioon on selline, mille puhul on sõltuval tunnusel ainult 2 taset (nt kas inimesel esineb konkreetne haigus või mitte). Multinomiaalse logistilise regressiooni puhul on sõltuva muutuja tasemeid rohkem kui 2.

Binaarne logistiline regressioon - ühe prediktoriga mudel

Andmetabel nimega “cowles” käsitleb seost tudengite isiksuseomaduste ja psühholoogilistes uurimustes osalemise valmiduse vahel. Andmestik koosneb neljast muutujast:

• neuroticism - Eysencki isiksuseküsimustku neurootilisuse alaskaala skoor;
• extraversion - ekstravertsuse alaskaala skoor;
• sex - sugu (female, male);
• volunteer - kas tudeng on valmis osalema edasises uurimistöös (no, yes);

Esmase ülevaate saamiseks andmetest kasutage funktsiooni summary.

summary(cowles)

##   neuroticism     extraversion       sex      volunteer
##  Min.   : 0.00   Min.   : 2.00   female:780   no :824  
##  1st Qu.: 8.00   1st Qu.:10.00   male  :641   yes:597  
##  Median :11.00   Median :13.00                         
##  Mean   :11.47   Mean   :12.37                         
##  3rd Qu.:15.00   3rd Qu.:15.00                         
##  Max.   :24.00   Max.   :23.00

Teeme kõigepealt logistilise regressioonmudeli, milles sõltuvaks tunnuseks on valmidus osaleda uurimistöös (tunnus volunteer) ja ennustame seda ekstravertsuse skoori kaudu. Logistilise mudeli tegemiseks kasutame funktsiooni glm (generalized linear model). Selle funktsiooni kasutamine on sarnane eelmises praktikumis kasutatud lm funktsiooniga. Kõigepealt peame sisestama mudelisse kaasatavate muutujate nimed (“sõltuv tunnus ~ sõltumatu tunnus”) ning seejärel andmestiku nime (“data = andmestik”). Funktsiooni glm puhul lisandub argument family, mille abil määrame, millist tüüpi mudelit teha tahame. Logistilise mudeli puhul paneme argumendi family väärtuseks binomial(). (Kui meil juhtub esinema puuduvaid andmeid, saame need välja jätta lisades argumendi na.action=na.omit. Samuti saaksime argumendi subset abil määrata valimit kitsendava tingimuse.)

cowles.mudel1 <- glm(volunteer ~ extraversion, data=cowles, family=binomial())

Uurime mudeli väljundit.

summary(cowles.mudel1)

## 
## Call:
## glm(formula = volunteer ~ extraversion, family = binomial(), 
##     data = cowles)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3379  -1.0584  -0.9299   1.2725   1.6243  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.13942    0.18538  -6.146 7.93e-10 ***
## extraversion  0.06561    0.01414   4.640 3.49e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1933.5  on 1420  degrees of freedom
## Residual deviance: 1911.5  on 1419  degrees of freedom
## AIC: 1915.5
## 
## Number of Fisher Scoring iterations: 4

Vaatame väljundi alaosa nimega Coefficients. Tulbas Estimate on ära toodud mudeli parameetrite - vabaliikme (Intercept) ja regressioonikordaja - väärtused. Tulbas Std. Error on ära toodud standardviga (kui suurt kordaja kõikumist võib oodata erinevates valimites). Tulbas z value on toodud z-statistiku väärtus, mis on tuntud ka Waldi statistiku nime all. See on saadud jagades regressioonikordaja väärtuse ja standardvea omaga. Tulbas Pr(>|z|) on märgitud z-statistiku kaudu arvutatud p-väärtus, mille abil saame hinnata, kas prediktor omab olulist seost sõltuva muutujaga. Praegusel juhul on ekstravertsuse regressioonikordajale vastav p-väärtus \(3.49e-06\) ehk \(3.49 * 10^{-6}\). Antud prediktori olulisusele nivool p < .001 viitavad ka rea lõpus olevad 3 tärni.

exp(coef(cowles.mudel1))

##  (Intercept) extraversion 
##     0.320005     1.067813

exp(confint(cowles.mudel1))

## Waiting for profiling to be done...

##                  2.5 %    97.5 %
## (Intercept)  0.2217862 0.4589052
## extraversion 1.0387831 1.0980308

cowles.nullmudel <- glm(volunteer ~ 1, data=cowles, family=binomial())

anova(cowles.nullmudel, cowles.mudel1, test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: volunteer ~ 1
## Model 2: volunteer ~ extraversion
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      1420     1933.5                          
## 2      1419     1911.5  1   22.022 2.695e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

install.packages("fmsb")

library(fmsb)
NagelkerkeR2(cowles.mudel1)

## $N
## [1] 1421
## 
## $R2
## [1] 0.02068324

cowles.mudel1 <- glm(volunteer ~ scale(extraversion), data=cowles, family=binomial())
exp(coef(cowles.mudel1))

##         (Intercept) scale(extraversion) 
##           0.7206583           1.2910767

Binaarne logistiline regressioon - kahe prediktoriga mudel

Mudelitesse saab prediktorina kaasata ka kategoriaalseid muutujaid. Teeme uue mudeli, lisades täiendava pretiktorina vastaja soo (tunnus sex) ja uurime mudeli väljundit.

cowles.mudel2 <- glm(volunteer ~ scale(extraversion) + sex, data=cowles, family=binomial())
summary(cowles.mudel2)

## 
## Call:
## glm(formula = volunteer ~ scale(extraversion) + sex, family = binomial(), 
##     data = cowles)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3860  -1.0499  -0.9035   1.2533   1.6853  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -0.21765    0.07259  -2.998  0.00271 ** 
## scale(extraversion)  0.25494    0.05519   4.619 3.85e-06 ***
## sexmale             -0.24662    0.10929  -2.256  0.02404 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1933.5  on 1420  degrees of freedom
## Residual deviance: 1906.4  on 1418  degrees of freedom
## AIC: 1912.4
## 
## Number of Fisher Scoring iterations: 4

Tabelist Coefficients näeme p-väärtuste abil, et mõlemad prediktorid on olulised. Kategoriaalsete muutujate puhul valitakse üks kategooriatest baaskategooriaks, millega ülejäänud kategooriaid võrreldakse. Praegusel juhul on muutujal sex kaks taset (female ja male). Vaikimisi on baaskategooriaks valitud female, tulenevalt sellest, et female on tähestukuliselt eespool võrreldes male’iga. Seega näitab vastav regressioonikordaja, kui erinev on kategooria male võrreldes kategooriaga female. Miinusmärgist selle kordaja ees võime järeldada, et meeste puhul on jaatava vastuse tõenäosus naistega võrreldes madalam. Kordajate arusaadavamaks tõlgendamiseks arvutame neist jällegi riskisuhted.

exp(coef(cowles.mudel2))

##         (Intercept) scale(extraversion)             sexmale 
##           0.8044055           1.2903863           0.7814391

Ekstravertsuse riskisuhe on jäänud samaks. Vastaja soo riskisuhtest võime järeldada, et meeste puhul on jaatava vastamise tõenäosus 78% naiste jaatava vastuse tõenäosusest. Ehk naiste puhul jaatava vastuse tõenäosus \(1/0.78 = 1.28\) korda suurem kui meeste puhul.

Hindame prediktoreid ka 95%-usaldusvahemike abil.

exp(confint(cowles.mudel2))

## Waiting for profiling to be done...

##                         2.5 %    97.5 %
## (Intercept)         0.6974271 0.9271049
## scale(extraversion) 1.1587831 1.4388556
## sexmale             0.6304839 0.9678052

Nagu näha ei sisalda kummagi prediktori vahemikud väärtust 1 ja sellest tulenevalt võib neid lugeda usaldusväärseteks.

Arvutame mudeli pseudo-determinatsioonikordaja.

NagelkerkeR2(cowles.mudel2)

## $N
## [1] 1421
## 
## $R2
## [1] 0.02543564

See on jäänud enam-vähem samaks, nii et soo lisamine mudeli kirjeldusvõimet väga palju paremaks ei teinud.

Multinomiaalne logistiline regressioon

Vaatame põgusalt ka logistilist regressiooni sõltuva muutuja korral, millel on rohkem kui kaks taset. Kasutame selleks paketi nnet funktsiooni multinom().

install.packages("nnet") # installeerime paketi

library(nnet) 

Andmetabelis ESS on 4 muutujat 2008. aasta Euroopa Sotsiaaluuringu Eesti vastajaid puudutavatest andmetest. üks tunnustest puudutab seda, millise erakonna poolt hääletas vastaja viimastel valimistel (tunnus partei). Näitlikustamise lihtsuse huvides olen antud tabelisse alles jätnud ainult 3 valimis kõige suuremat toetust omavat parteid (IRL, Kesk ja Reform). Teeme mudeli, milles ennustame erakondlikku eelistust vastaja vanuse kaudu. Vaikimisi võetakse baaskategooriaks tähestiku järjekorras kõige esimene kategooria, antud juhul oleks selleks IRL. (Baaskategooria saab muuta funktsiooni relevel abil.)

ess.mudel1 <- multinom(partei ~ vanus, data = ESS)

## # weights:  9 (4 variable)
## initial  value 727.281335 
## final  value 677.750378 
## converged

Uurime tulemusi:

summary(ess.mudel1)

## Call:
## multinom(formula = partei ~ vanus, data = ESS)
## 
## Coefficients:
##        (Intercept)        vanus
## Kesk    -0.6269266  0.023066461
## Reform   1.1064595 -0.006418908
## 
## Std. Errors:
##        (Intercept)       vanus
## Kesk     0.3475546 0.006384485
## Reform   0.3145498 0.006113478
## 
## Residual Deviance: 1355.501 
## AIC: 1363.501

Väljund on jaotatud kaheks osaks: regressioonikordajad ja standardvead. Regressioonikordajad on toodud kummagi partei jaoks eraldi. Nende põhjal on näha, et võrreldes IRL-iga suurendab kasvav vanus Keskerakonna poolt hääletamise tõenäosust, samas kui Reformierakonna toetamise tõenäosusega ei tundu seost olevat. Kahjuks ei anna multinom funktsioon meile p-väärtusi. Nende välja arvutamiseks on erinevaid viise, aga kõige lihtsam on kasutada paketi stargazer samanimelist funktsiooni, mis annab meile regressioonimudeli tulemused ja lisab sinna p-väärtused.

install.packages("stargazer")

library(stargazer)
stargazer(ess.mudel1, type="text", out="multi1.htm")

## 
## ==============================================
##                       Dependent variable:     
##                   ----------------------------
##                        Kesk         Reform    
##                        (1)            (2)     
## ----------------------------------------------
## vanus                0.023***       -0.006    
##                      (0.006)        (0.006)   
##                                               
## Constant             -0.627*       1.106***   
##                      (0.348)        (0.315)   
##                                               
## ----------------------------------------------
## Akaike Inf. Crit.   1,363.501      1,363.501  
## ==============================================
## Note:              *p<0.1; **p<0.05; ***p<0.01

Vanuse mõju suurust saab hinnata riskisuhete abil:

exp(coef(ess.mudel1))

##        (Intercept)     vanus
## Kesk     0.5342312 1.0233345
## Reform   3.0236343 0.9936016

Keskerakonna regressioonikordaja riskisuhte alusel võime järeldada, et lisanduv eluaasta suurendab Keskerakonna toetamise suhtelist tõenäosust võrreldes IRL-i toetamise tõenäosusega 1.02 korda ehk 2%.

Ülesanded

1. Andmetabelist nimega neeme on Marko Neeme magistritöö (Neeme, 2012) andmed, milles uuriti 50-70-aastaste meeste suhtumist eesnäärmevähi skriiningtesti. Tunnus “”valmidus“, näitab testis osalemise valmidust (tasemed pigem jah, ja pigem ei). Lisaks on tabelis ära toodud Suure Viisiku isiksuseomadused: neurootilisus (tunnus N), ekstravertsus (E), avatus kogemusele (O), sotsiaalsus (A) ja meelekindlus (C). Koostage binaarse logistilise regressiooni mudel, mis ennustab skriiningtestis osalemise valmidust Suure Viisiku isiksuseomaduste kaudu.

1. Esialgu tuleks saada ülevaade andmetest. Kasutage funktsiooni summary. Milline on oslejate keskmine vanus? Missugune on “valmiduse” jaotus?

2. Tehke ekstravertsuse tulemustest histogramm.

3. Millised isiksusomadused omavad olulist seost valmidusega? Kas need omadused suurendavad või vähendavad testis osalemise valmidust?

4. Arvutage välja riskisuhted ja nende 95%-usaldusvahemikud? Milline omadus mõjutab osalemisvalmidust kõige tugevamini?

5. Arvutage mudeli sobitusastet näitav hii-ruut-statistik ja selle p-väärtus.

6. Leidke mudeli pseudo-determinatsioonikordja.

LISAD

cdplot(volunteer ~ extraversion, data=cowles)

expit <- function(x) exp(x) / (1+exp(x))

coef(cowles.mudel1)

##         (Intercept) scale(extraversion) 
##          -0.3275901           0.2554766

expit(-0.3275901 + 2*0.2554766)

## [1] 0.5457128