R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Ordinaalinen ja multinomiaalinen regressio

Jos vastemuuttuja on kaksiluokkainen (esim. naaras/koiras), käytetään aineiston analysointiin tyypillisesti logistista regressiota. Ordinaalinen logistinen regressio on logistisen regression laajennus siinä mielessä että se on menetelmä, jota voidaan käyttää, kun vastemuuttuja on moniluokkainen ja lisäksi järjestysasteikollinen (esim. vähän/kohtalaisesti/paljon). Keskeinen oletus on ”proportional odds” tai ”rinnakkaisten regressioiden oletus” eli oletetaan että vastemuuttujan tasojen muutos vaikuttaa samalla tavalla kaikkiin sovitettuihin logit-funktioihin. Jos oletus ei päde, voidaan käyttää multinominaalista regressiota, joka ei tee samaa oletusta, mutta jättää myös vastemuuttujan tasojen järjestyksen huomiotta.

Käytetäänpä menetelmien havainnollistamiseksi esimerkkinä European Values Studyn (EVS) dataa Suomesta viimeisimmältä kierrokselta (2008-2010, Suomesta vuodelta 2009). ESV:n aineisto on saatavissa vain akateemista tutkimusta varten. Tämä kirjoitus täyttää käyttöehdon, sillä se on syntynyt tutkimusprojektissa osana aineiston katselmointia, muttei kata koko tutkimusta eikä ole lopullinen julkaisu.

Aineiston lataaminen R:ään

Aineiston lataaminen käy helpoimmin SPSS-tiedostosta, sillä foreign-paketti, joka tarjoaa sopivan funktion, tulee oletusarvoisesti R-asennuksen mukana. Toinen vaihtoehto on käyttää pakettia sas7bdat, joka lukee SAS-muotoisia tiedostoja. Aineisto, jossa on vain Suomen kyselyn tulokset vuodelta 209, saadaan luettua SPSS_tiedostosta seuraavasti:

setwd("C:\\Users\\tuimala_j\\Desktop\\EVS\\fi2008")
library(foreign)
d<-read.spss("ZA4762_v1-1-0.sav", to.data.frame=TRUE)

Tuloksena on data frame, jossa on 1134 riviä ja 448 saraketta.

Aineiston muokkaus

Aineistoa muokattiin enemmän analyyseihin sopivaan muotoon. Esimerkiksi, henkilöille laskettiin ikä, ja monien muuttujien koodausta muutettiin siten, että pienimmät vastaajaryhmät yhdistettiin isompiin, muuttujien suunta tietysti huomioiden. Alle on koottu pääasialliset muokkaukset:

d$age<-2009-d$v303
d$gender<-d$v302
d$taxcheat<-d$v234 # 1=never, 10=always
d$taxgrey<-d$v245  # 1=never, 10=always
d$wouldvote<-d$v263
d$maritalstatus<-as.character(d$v313)
d$maritalstatus[as.numeric(d$v313)<=2]<-"married"
d$maritalstatus[as.numeric(d$v313)>=3 & as.numeric(d$v313)<=5]<-"not married"
d$maritalstatus[as.numeric(d$v313)==6]<-"never married"
d$maritalstatus<-as.factor(d$maritalstatus)
d$education<-as.character(d$v336)
d$education[as.numeric(d$v336)<=3]<-"low"
d$education[as.numeric(d$v336)>=4 & as.numeric(d$v336)<=5]<-"middle"
d$education[as.numeric(d$v336)>=6]<-"high"
d$education<-as.factor(d$education)
d$trustparliament<-d$v211
levels(d$trustparliament)<-c("quite a lot", "quite a lot", "not very much", "none at all")
d$trustgovernment<-d$v222
levels(d$trustgovernment)<-c("quite a lot", "quite a lot", "not very much", "none at all")
d$trustcivilservants<-d$v212
levels(d$trustcivilservants)<-c("quite a lot", "quite a lot", "not very much", "none at all")
d$satdemoc<-d$v223
levels(d$satdemoc)<-c("satisfied", "satisfied", "not satisfied", "not satisfied")
d$proud<-d$v256
levels(d$proud)<-c("very proud", "quite proud", "not proud", "not proud")
d$happy<-d$v8
levels(d$happy)<-c("very happy", "quite happy", "not happy", "not happy")
d$health<-d$v9
levels(d$health)<-c("good health", "good health", "fair health", "poor health", "poor health")
d$trustpeople<-d$v62
d$lifesatisfaction<-d$v66
d$politicalleaning<-d$v193
d$incomecat<-d$v353YR
levels(d$incomecat)<-c("poor", "poor", "poor", "poor", "low", "low", "low", "middle", "middle", "high", "high", "high")
d$employment<-d$v337
levels(d$employment)<-c("employed", "employed", "self employed", "not employed", "pensioned", "not employed", "not employed", "not employed", "not employed", "not employed")
d$reforms<-d$v200
d$religiousness<-d$v114
d$morals<-d$v104

Puuttuvien havaintojen imputointi

Analyysejä varten aineistosta erotettiin osa-aineisto, jonka puuttuvat havainnot imputoitiin mice-paketilla. Saatu imputoitu aineisto ei ole täysin optimaalinen, ja vaatisi vielä lisätyötä, jotta siitä saisi paremman. Tällä tavoin imputoitu aineisto antaa kuitenkin melko tarkkaan samoja tuloksia kuin imputoimatonkin, js riittänee siten ainakin tähän esimerkkiin. Imputointi tehtiin seuraavasti:

d2<-d[,colnames(d) %in% c("taxcheat", "age","gender","maritalstatus","education","employment","incomecat","wouldvote","reforms","religiousness","morals","happy","proud","health","trustparliament","trustgovernment","trustcivilservants","weight")]
library(mice)
d2.1<-d2
d2.1$taxcheat<-factor(d2.1$taxcheat, ordered=TRUE)
d3<-complete(mice(d2.1, method="polr", m=1))
d4<-na.omit(d2)

Havaitun aineiston ja imputoidun aineiston vertaaminen imputoinnin onnistumisen selvittämiseksi onnistuu esimerkiksi graafisesti tai vaihtoehtoisesti taulukoimalla. Muuttujassa reforms, joka mittaa perusideologista asennetta, on kaikkein eniten puuttuvia havaintoja (124 kpl). Verrataanpa sen jakaumaa aineistossa, josta on poistettu puuttuvat havainnot (d4, punaisella) ja imputoidussa aineistossa (d3, mustalla):

# Muuttuja on aineiston sarakkeessa 16
# Tiheysfunktio
plot(density(as.numeric(d3[,16])))
lines(density(as.numeric(d4[,16])), col="red")

# Taulukko
# Prosenttiosuuksien ero aineistojen välillä
(prop.table(table(d4[,16]))-prop.table(table(d3[,16])))*100
 
radically changed by revolutionary action              gradually changed by reforms 
                                0.2204586                                -2.0655271 
             defended against all changes 
                                1.8450685

Sekä taulukosta että kuvasta on mahdollista päätellä, että alkuperäisen ja imputoidun aineiston jakauma reform-muuttujan osalta on varsin samanlainen, eroa frekvensseissä on korkeintaan 2%. Samalla tavalla voidaan havaita, että useimmat muuttujat ovat samaan tapaa varsin samanlailla jakautuneita sekä alkuperäisessä että imputoidussa aineistossa, mikä on hyvä.

Graafisia tarkasteluja

Aineistossa on useita, erilaisten asioiden hyväksyttyvyyttä luotaavia kysymyksiä, joista monet ovat potentiaalisesti kiinnostavia vastemuuttujia. Vastaukset ovat kaikissa asteikolla 1-10, jossa 1 merkitsee ”Ei koskaan hyväksyttävää” ja 10 ”Aina hyväksyttävää”. Muuttujia voidaan tarkastella graafisesti yhdessä esimerkiksi seuraavasti. Ensin muuttujista muodostetaan uusi likert-tyyppinen aineisto, jonka muodostaminen onnistuu vain, jos kaikissa muuttujissa on sama mitta-asteikko. Tämän jälkeen muuttujat voidaan visualisoida pylväskaaviona:

library(likert)
test<-na.omit(d[,250:269])
colnames(test)<-paste(colnames(test), c("claim state benefits",
                  "cheating on tax",
                  "joyriding",
                  "taking soft drugs",
                  "lying in own interest",
                  "adultery",
                  "accepting a bribe",
                  "homosexuality",
                  "abortion",
                  "divorce",
                  "euthanasia",
                  "suicide",
                  "paying in cash to avoid taxes",
                  "having casual sex",
                  "avoiding public transportaion fare",
                  "prostitution",
                  "experiments using human embryos",
                  "mainpulation of food",
                  "in vitro fertilization",
                  "death penalty"))
li<-likert(test)
library(HH)
plot(li)
detach(package:HH)

Muodostuva kuvio on seuraavankaltainen:

Suurin osa vastaajista tuomitsee esimerkiksi verovilpin, sosiaalietuuksien huijaamisen ja julkisessa liikenteessä matkustamisen ilman lippua. Seuraavissa analyyseissä käytetään vastemuuttujana verovilppiä (v234).

Muiden muuttujien suhdetta vastemuuttujiin voidaan tarkatella useallakin tavalla, mutta eräs varsin nopeasti silmäiltävä on tavanomainen pinottu pylväskaavio. Esimerkiksi sukupuolen vaikutusta veronmaksumyönteisyysteen voidaan tarkastella pinoamalla eri mielipiteet sukupuolittain pylväskaavioon:

library(RColorBrewer)
cols<-c(brewer.pal(9, "Reds"), "black")
barplot(prop.table(table(d4[,4], d4[,3]), 2), col=(cols), border=(cols), axes=FALSE, main="Cheating on tax")
abline(h=seq(0.1,1,0.1), col="white")
axis(side=2, at=seq(0,1,0.1), labels=seq(0,1,0.1), las=1)

Tuloksena on seuraavanlainen pylväskaavio, josta on melko nopeasti nähtävissä, että naiset ovat veronmaksumyönteisempiä kuin miehet. Miehistä n. 85% on voimakkaasti tai melko voimakkaasti sitä mieltä, että huijaaminen veroissa ei ole hyväksyttävää, mutta naisista tätä mieltä on vastaavasti jopa n. 95%. Edellinen kuva on laadittu alkuperäisestä aineistosta, mutta imputoidulla aineistollakin kuva on oleellisesti samanlainen.

Samaan tapaan voidaan havaita esimerkiksi luottamuksen virkamiehiin vaikuttavan veronmaksumyönteisyyteen.

Tilastollinen analyysi

Painotus

Koska kyseessä on survey, jossa on käytetty suunniteltua otantaa, pitäisi analyyseissä käyttää myös painotusta. Paino on aineistossa mukana nimellä weight. Kaikkein vanhimmissa ikäryhmissä painot alkavat olla melko suuria, sillä yli 76-vuotiaitten ryhmässä kukin kyselyyn osallistunut vastaa karkeasti n. 15 000 – 30 000 asukasta, mikä voi vaikuttaa tulosten luotettavuuteen. Jos oiotaan hieman mutkia, ja määritellään otanta-asetelma R:n survey-paketin avulla seuraavasti, nähdään pääperiaate painojen käyttämisessä:

library(survey)
d3$strat<-factor(d3$weight, labels=letters[1:14])
d3$weight2<-d3$weight*5351427/1134
s<-svydesign(id=~1, weights=~weight2, strata=~strat, data=d3)

Yllä on muutettu analyysipainot (weight) asetelmapainoiksi (weight2) totaalien estimointia varten, ja tehty oletus että otanta on stratifioitu iän ja sukupuolen mukaan, joka surveyn dokumentaation mukaan ei välttämättä ole aivan tarkka lähestymistapa, mutta EVS on laskenut painot näissä (jälki)ositteissa (?). Näitä käytetään sitten otanta-asetelman kuvaamiseksi svydesign()-funktiolla. Näin muodostettua aineistoa käyttäen voidaan arvioida esimerkiksi eri verovilppiluokkiin kuuluvien asiakkaiden määrä koko populaatiossa:

svytotal(~taxcheat, s)
             total    SE
taxcheat1  3430581 83454
taxcheat2   799569 63809
taxcheat3   500452 52676
taxcheat4   193180 29741
taxcheat5   143723 24937
taxcheat6   113240 22727
taxcheat7    80677 19038
taxcheat8    36658 12203
taxcheat9    22662 10321
taxcheat10   30685 11741

Näin ollen näyttää siltä, että väestöstä n. 64% (3,4 miljoonaa asukasta) on sitä mieltä, että verovilppi ei ole milloinkaan perusteltavissa tai sallittua. Tässä painot on skaalattua koko populaatioon, vaikka surveyssä mukana olivat vain 18-74 vuotiaat asukkaat. Jos käytetään ko. ikäryhmiin perustuvaa skaalausta, on täysi veromyönteisiä väestössä sama prosenttiosuus, mutta lukumääräisesti n. 2,5 miljoonaa.

Edellä esitetty on parasta tulkita suuntaa-antavana, koska painotuksen käsittely ei ole tässä täysin oikein toteutettu. Osittain tämä johtuu siitä, että dokumentaatio on tässä kohden hieman epäselvää, eikä kaikki haluamiani painotuksen laskemiseen tarvittavia tietoa ole dokumenaatiossa. Niinpä lopulliseen analyysin painotuksen pohtiminen pitää tehdä vielä uudelleen. Lisäksi kaikki seuraavassa tarvittavia menetelmiä ei löydy survey-paketista, ja painojen käyttö analyyseissä vaatisi lisätyötä. Esimerkiksi multinomiaalisessa regressiossa varten pitäisi generoida uusi painotus.

Lineaarinen regressio

Analyseissä tarkasteltava vastemuuttuja voi saada kymmenen erillistä arvoa, joten voisi ajatella, että analyysin voisi tehdä lineaarista regressiota käyttäen. Idea on helppo testata:

fit<-(lm(taxcheat~age+gender+maritalstatus+education+reforms+religiousness+morals+happy+proud+health+trustparliament+trustgovernment+trustcivilservants, data=d3))
par(mfrow=c(2,2))
plot(fit)

Yläpuolella esitetty lineaarisen regression oletusten tarkistamiseen tarkoitettu kuvapatteri antaa selkeän kuvan siitä, että lineaarinen regressio ei liene menetelmä, joka tässä tapauksessa tulisi kyseeseen. Mallin residuaalinen jakauma on monessakin suhteessa niin patologinen, ettei sitä kannattane alkaa sen paremmin selvittelemään tai ongelmia ratkomaan.

Lineaarisen regression sijaan kyseeseen voi tulla ordinaalinen regressio.

Ordinaalinen regressio

Ordinaalisen regression sovittamiseen soveltuvia paketteja ja funktioita on tarjolla useita. Tässä käytetään ordinal-paketin tarjoamia mahdollisuuksia. Malli voidaan sovittaa seuraavasti:

library(ordinal)
or<-clm(taxcheat~age+gender+maritalstatus+education+reforms+religiousness+morals+happy+proud+health+trustparliament+trustgovernment+trustcivilservants, data=d8)

Ennen mallin sen tarkempaa tarkastelua on kuitenkin syytä testata mallin oletuksia. Tämä tapahtuu vertaamalla kahta erilaista mallia, yksi selittävä muuttuja kerrallaan. Sovitetaan siis kaksi mallia vaikkapa sukupuoli-muuttujan oletusten tarkastamiseksi:

or1<-clm(taxcheat~age+gender, data=d8)
or2<-clm(taxcheat~age, nominal=~gender, data=d8)
anova(or1, or2)
 
Likelihood ratio tests of cumulative link models:
 
    formula:                nominal: link: threshold:
or1 taxcheat ~ age + gender ~1       logit flexible  
or2 taxcheat ~ age          ~gender  logit flexible  
 
    no.par    AIC  logLik LR.stat df Pr(>Chisq)   
or1     11 2859.6 -1418.8                         
or2     19 2851.4 -1406.7  24.208  8   0.002115 **
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Koska testin tulos (P=0.002115) on merkitsevä, ei proportional odds-oletus päde sukupuoli-muuttujalle. Koska aikaisemmin olemme todenneet, että sukupuoli oleellinen huomioitava muuttuja, ei myöskään ordinal regression taida olla tälle aineistolle sopiva menetelmä tai lähestymistapa. Joudumme siis turvautumaan ns. viimeiseen oljenkorteen multinomiaaliseen regressioon.

Multinomiaalinen regressio

Multinomiaalinen regressio sovitetaan nnet-paketin funktiolla multinom(). Voimme aluksi sovittaa yksinkertaisen malli, jossa ovat vain ikä ja sukupuoli:

mt<-multinom(taxcheat~age+gender, data=d3)
summary(mt)
Call:
multinom(formula = taxcheat ~ age + gender, data = d3)
 
Coefficients:
   (Intercept)         age     gender2
2   -0.7209358 -0.01476196  0.01911904
3   -0.2045480 -0.03218150 -0.42059019
4   -1.6115846 -0.02103312 -0.48683715
5   -1.1209869 -0.03422574 -0.91428835
6   -0.8345461 -0.05439850 -0.33137491
7   -1.2433190 -0.04046795 -1.84105549
8   -2.2221593 -0.03525719 -1.47590966
9   -2.3393510 -0.04604764 -1.62549101
10  -1.8291585 -0.04480941 -8.36080279
 
Std. Errors:
   (Intercept)         age      gender2
2    0.2978832 0.005705915 0.1707849160
3    0.3443524 0.007143866 0.2091671250
4    0.5337961 0.010755660 0.3218782910
5    0.5730468 0.012344323 0.3787640435
6    0.6490525 0.015061843 0.4107449074
7    0.7420633 0.016741331 0.6386489196
8    1.0517747 0.023167003 0.8068195621
9    1.3696389 0.031834783 1.1220843009
10   1.1445682 0.026748266 0.0004675733
 
Residual Deviance: 2803.023 
AIC: 2857.023

Kullekin vastemuuttujan tasolle on sovitettu oma funktio, ja selittäville tekijöille estimoidut vaikutusten arviot voidaan tulkita log odds-asteikolla, kuten tavanomaisen logistisen regression arviotkin. Siten relative risk kuulua vastemuuttujan luokkaan 10 on exp(-8.36080279)~0.0002 verrattuna luokkaan 1, jos enkilö on nainen. Toisin sanoen, veromyönteisyys on naisilla miehiä suurempaa (koska pieni vastemuuttujan vastasi suurempaa myöntyvyyttä). Voimme seuraavaksi sovittaa monimutkaisemman mallin esimerkiksi seuraavasti:

mt2<-multinom(taxcheat~age+gender+maritalstatus+education+reforms+religiousness+morals+happy+proud+health+trustparliament+trustgovernment+trustcivilservants, data=d3)

Tuloksena on kuitenkin tässä vaikeasti tulkittava malli, koska vastemuuttujassa on niin monia tasoja, ja osassa tasoista vain hyvin harvoja havaintoja. Malleja voidaan verrata vaikkapa AIC-arvoilla (funtio AIC()), jolloin havaintaan, että monimutkaisempi malli on yksinkertaisempaa huonommin aineistoon sopiva.

Luodaanpa aineistoon ensin uusi, kolmiluokkainen muuttuja, joka kuvastaa veromyönteisyyttä, ja sovitetaan sitten aiemmin esitetyt mallit uudelleen:

d3$taxcut<-cut(as.numeric(d3$taxcheat), 3)
mt<-multinom(taxcut~age+gender, data=d3)
mt2<-multinom(taxcut~age+gender+maritalstatus+education+reforms+religiousness+morals+happy+proud+health+trustparliament+trustgovernment+trustcivilservants, data=d3)
AIC(mt)
[1] 731.0067
AIC(mt2)
[1] 751.8448

Monimutkaismepi malli on edelleen huonommin sopiva, mutta nyt tulokset ovat jo helpommin tulkittavia. Muuttujien vaikutuksille voidaan laskea myös normaaliapproksimaatiolla p-arvot, jos ne on aivan pakko laskea:

# t-arvot
z <- summary(mt)$coefficients/summary(mt)$standard.errors
# p-arvot
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p
       (Intercept)          age     gender2
(4,7]  0.050623362 3.557384e-05 0.001175424
(7,10] 0.007404069 2.982090e-02 0.001859389

Kuten tuloksesta huomataan ovat iän ja sukupuolen vaikutukset merkitseviä. Sovitetaanpa vielä lopuksi yksi aiempia paremmin sopiva malli, ja tehdään sille muutamia testejä, kuten edelläkin:

mt<-multinom(taxcut~age+gender+trustcivilservants+proud+happy, data=d3)
summary(mt)
 
#Call:
#multinom(formula = taxcut ~ age + gender + trustcivilservants + 
#    proud + happy, data = d3)
#
#Coefficients:
#       (Intercept)         age    gender2 trustcivilservants2 trustcivilservants3     proud2    proud3
#(4,7]    -1.822202 -0.03482159 -0.7676242           0.5720857          -0.2679447  0.1822489 0.5871784
#(7,10]   -2.714720 -0.03504897 -1.8915358           1.3410950           2.1293038 -0.6617611 0.6747209
#           happy2    happy3
#(4,7]   0.6299150 1.2840112
#(7,10] -0.1267706 0.2342456
#
#Std. Errors:
#       (Intercept)        age   gender2 trustcivilservants2 trustcivilservants3    proud2    proud3
#(4,7]    0.5775368 0.00859132 0.2551511           0.2650046           0.4787611 0.2620872 0.4243382
#(7,10]   1.0414296 0.01670200 0.6314767           0.6611523           0.7169703 0.5560131 0.6492743
#          happy2    happy3
#(4,7]  0.4488657 0.5047430
#(7,10] 0.6775498 0.8206067
#
#Residual Deviance: 684.567 
#AIC: 720.567 
 
z <- summary(mt)$coefficients/summary(mt)$standard.errors
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p
 
#       (Intercept)          age     gender2 trustcivilservants2 trustcivilservants3    proud2    proud3
#(4,7]  0.001604283 5.054072e-05 0.002625338          0.03086748         0.575709581 0.4868203 0.1664348
#(7,10] 0.009141288 3.586193e-02 0.002740704          0.04251733         0.002979325 0.2339718 0.2987153
#         happy2     happy3
#(4,7]  0.160513 0.01096253
#(7,10] 0.851581 0.77529621
 
exp(summary(mt)$coefficients)
 
#       (Intercept)       age   gender2 trustcivilservants2 trustcivilservants3    proud2   proud3
#(4,7]   0.16166939 0.9657777 0.4641144            1.771959           0.7649501 1.1999128 1.798905
#(7,10]  0.06622349 0.9655581 0.1508400            3.823228           8.4090101 0.5159419 1.963485
#          happy2   happy3
#(4,7]  1.8774511 3.611095
#(7,10] 0.8809357 1.263955

Mallin perusteella tulkittuna ihmiset, jotka eivät luota virkamiehiin juuri lainkaan (trustcivilservants3) ovat myös kaikkein vähiten veronmaksumyönteisiä. Lisäksi mm. ikä ja sukupuoli ovat myös ennustavia tekijöitä, mutta paljon heikompia kuin luottamus virkamiehiin. Mielenkiintoisesti ylpeys omasta maasta (proud) on melko voimakas ennustava tekijä. Vähäisempi ylpeys on assosioitunut vähäisempään veromyönteisyyteen.

Lisätietoa menetelmistä

Seuraavissa dokumenteissa on esitetty eri malleja ja niiden sovittamista ja tarkastelua hieman ylläolevaa tarkemmalla tasolla:


Vastaa

Sähköpostiosoitettasi ei julkaista. Pakolliset kentät on merkitty *