Fysisk samling 05.03. Multippel og logistisk regresjon

Njål Foldnes

Multippel og logistisk regresjon

Multippel regresjon er et hovedverktøy i psykologi

  • Fra forskningsspørsmålet identifiserer vi en avhengig variabel \(Y\) som vi vil forklare variasjonen til
  • Og en rekke prediktorer \(X_1, X_2, \ldots, X_k\) som skal forklare variasjonen i \(Y\)
  • Modellen:

\[ Y = \beta_0 + \beta_1 X_1+ \beta_1 X_2+\ldots + \beta_k X_k+\epsilon\]

Fra forskningsspørsmål til regresjon

  • Finne de viktigste prediktorene
  • Sjekke forutsetninger: stort nok utvalg, unngå multikollinearitet1, sammenhengene er lineære, homoskedastisitet2, normalfordelte residualer
  • Sjekke modellens totale forklaringskraft (Adjusted \(R^2\) og \(F\))
  • Sjekke hver prediktors forklaringskraft
  • Besvare forskningsspørsmålet

Digitalt stress

Vi ser på noen data fra stress-studien. Dataene er i comma separated format og lastes ned her

For å lese inn data i mange format (excel, csv, spss, osv) trykk på File>Import dataset i Rstudio

library("tidyverse")
stressdata <- read.csv("../data/stress_complete.csv") 
stressdata$X <- NULL#fjerne unødvendig variabel "X"
dim(stressdata)#rader og søyler
[1] 581  23
colnames(stressdata)#hvilke variabler
 [1] "age"                     "sex"                    
 [3] "civil_status"            "education"              
 [5] "professional_experience" "COMP1"                  
 [7] "COMP2"                   "COMP3"                  
 [9] "COMP4"                   "COMP5"                  
[11] "CONF1"                   "CONF2"                  
[13] "CONF3"                   "CONF4"                  
[15] "CONF5"                   "INSE1"                  
[17] "INSE2"                   "INSE3"                  
[19] "INSE4"                   "INSE5"                  
[21] "alene"                   "leder"                  
[23] "bransje"

Beskrivende statistikk - Alder

Vi begynner alltid med å se på grafer

hist(stressdata$age)

Vi har outlier. Gjennomsnitt 52.6 og median: 42. Fjern outlier:

stressdata <- filter(stressdata, age < 70)
nrow(stressdata)
[1] 576
hist(stressdata$age)

Beskrivende statistikk - Utdanningsnivå

Fire utdanningsnivå (grunnskole, vid.g, bachelor, master/phd)

qplot(stressdata$education)+geom_bar()

Andeler

prop.table(table(stressdata$education))

         1          2          3          4 
0.01388889 0.30034722 0.39409722 0.29166667

Beskrivende statistikk kjønn

Variabelen “sex” er kodet slik: 1=kvinne, 2=mann, 3=ikke oppgitt.

Vi forenkler med å lage en ny binær variabel “kjønn”

stressdata$kjønn <- ifelse(stressdata$sex == 2,"mann", "kvinne")
qplot(stressdata$kjønn)+geom_bar()

Beskrivende statistikk bransje

table(stressdata$bransje )

                    annet            Bygg og anlegg         Faglige tjenester 
                       47                        14                         5 
     Finans og forsikring Forskning og undervisning           Helse og omsorg 
                       10                        47                       347 
                 Industri                        IT  Offentlig administrasjon 
                        7                        23                        25 
             Olje og gass                Varehandel 
                       34                        17

Vi forenkler ved å lage ny binær variabel for bransje (vi overskriver den gamle bransje)

stressdata$bransje <- ifelse(stressdata$bransje =="Helse og omsorg", "Helse og omsorg", "Annet")
table(stressdata$bransje)

          Annet Helse og omsorg 
            229             347

Complexity ICT er et aspekt med digitalt stress

Fem items som tilhører “Complexity ICT” er bra korrelert

comp <- select(stressdata, contains("COMP"))#velge ut bare disse 5
cor(comp) %>% round(1)
      COMP1 COMP2 COMP3 COMP4 COMP5
COMP1   1.0   0.8   0.8   0.7   0.7
COMP2   0.8   1.0   0.8   0.7   0.7
COMP3   0.8   0.8   1.0   0.7   0.7
COMP4   0.7   0.7   0.7   1.0   0.8
COMP5   0.7   0.7   0.7   0.8   1.0

Oppsummere complexity med en summeskår

Lage summeskår som et mål på Complexity stress

stressdata$complexity <- stressdata$COMP1+stressdata$COMP2+stressdata$COMP3+stressdata$COMP4+stressdata$COMP5
summary(stressdata$complexity)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   5.00   11.00   19.00   18.66   26.00   35.00 
hist(stressdata$complexity)

Kan vi forklare noe i variasjonen til complexity stress?

Kan kjønn og alder forklare noe av stresset?

mod <- lm(complexity ~ kjønn+age, data=stressdata) 
summary(mod)

Call:
lm(formula = complexity ~ kjønn + age, data = stressdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.1812  -7.1336  -0.1302   6.9751  18.6830 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 17.22973    1.44498  11.924  < 2e-16 ***
kjønnmann   -2.91213    0.88412  -3.294  0.00105 ** 
age          0.04760    0.03322   1.433  0.15244    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.452 on 573 degrees of freedom
Multiple R-squared:  0.02046,   Adjusted R-squared:  0.01704 
F-statistic: 5.983 on 2 and 573 DF,  p-value: 0.002681
  • p-verdien til F-testen er signifikant. Det betyr at kjønn og alder forklarer mer enn ingenting
  • \(R^2\) er imidlertid veldig lav

Endre modell

Ta bort alder og ta heller med bransje

mod <- lm(complexity ~ kjønn+bransje, data=stressdata) 
summary(mod)

Call:
lm(formula = complexity ~ kjønn + bransje, data = stressdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.3573  -6.3573   0.3255   6.3696  19.3255 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)             14.6745     0.6052  24.247   <2e-16 ***
kjønnmann               -0.1764     0.8628  -0.204    0.838    
bransjeHelse og omsorg   6.6828     0.7071   9.452   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.875 on 573 degrees of freedom
Multiple R-squared:  0.1495,    Adjusted R-squared:  0.1466 
F-statistic: 50.38 on 2 and 573 DF,  p-value: < 2.2e-16
  • p-verdien til F-testen er signifikant og \(R^2\) er mye høyere nå
  • Nå er ikke kjønn signifikant lengre.

Helsesektor er full av kvinner..

Det er helsesektoren som forklarer stress! Grunnen til at kjønn var signifikant før var at det er mange damer i den stressfulle helsesektoren

stressdata %>% group_by(bransje) %>% summarise(andel_menn=mean(kjønn=="mann"))
# A tibble: 2 × 2
  bransje         andel_menn
  <chr>                <dbl>
1 Annet               0.358 
2 Helse og omsorg     0.0980

Bransje og lederansvar

Ta bort kjønn og ta med lederansvar

mod <- lm(complexity ~ leder+bransje, data=stressdata) 
summary(mod)

Call:
lm(formula = complexity ~ leder + bransje, data = stressdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.3492  -6.3492   0.3539   6.3871  19.3539 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)             14.6461     0.5786  25.314   <2e-16 ***
lederTRUE               -0.1327     0.9649  -0.138    0.891    
bransjeHelse og omsorg   6.7031     0.6958   9.633   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.875 on 573 degrees of freedom
Multiple R-squared:  0.1495,    Adjusted R-squared:  0.1465 
F-statistic: 50.36 on 2 and 573 DF,  p-value: < 2.2e-16

Leder er ikke signifikant.

Utdanningsnivå og stress

Utdanningsnivå er på ordinal skala så vi gjør det om til factor:

mod <- lm(complexity ~ education+bransje, data=stressdata) 
summary(mod)

Call:
lm(formula = complexity ~ education + bransje, data = stressdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.6117  -6.3274   0.2403   6.5246  19.3824 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)             15.0441     1.3507  11.138   <2e-16 ***
education               -0.1422     0.4095  -0.347    0.729    
bransjeHelse og omsorg   6.7098     0.6727   9.975   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.875 on 573 degrees of freedom
Multiple R-squared:  0.1497,    Adjusted R-squared:  0.1467 
F-statistic: 50.42 on 2 and 573 DF,  p-value: < 2.2e-16

Utdanningsnivå er heller ikke signifikant prediktor

Parsimony = Enkelhetsprinsippet

Ironisk nok ender vi opp med en enkel regresjon! Det er bare bransje som forklarer stress i våre data (15% er forklart, 85% er ikke forklart)

mod <- lm(complexity ~ bransje, data=stressdata) 
summary(mod)

Call:
lm(formula = complexity ~ bransje, data = stressdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.3401  -6.3401   0.3886   6.3886  19.3886 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)             14.6114     0.5200   28.10   <2e-16 ***
bransjeHelse og omsorg   6.7287     0.6699   10.04   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.869 on 574 degrees of freedom
Multiple R-squared:  0.1495,    Adjusted R-squared:  0.148 
F-statistic: 100.9 on 1 and 574 DF,  p-value: < 2.2e-16

Når \(Y\) er binær: Logistisk regresjon

  • I multippel regresjon er \(Y\) kontinuerlig
    • For eksempel \(Y\)= well-being på en skala fra 0-10
    • Modellen gir forventet well-being gitt feks kjønn og alder
  • Hvis \(Y\) er binær så må vi bruke logistisk regresjon
    • For eksempel \(Y\)=depressiv (Ja/Nei)
    • Modellen gir da odds/sannsynlighet for at en person er depressiv, gitt feks kjønn og alder

Eksempel: Extramarital affairs

library(AER)#install first if needed
data("Affairs")#loader datasettet
summary(Affairs)
    affairs          gender         age         yearsmarried    children 
 Min.   : 0.000   female:315   Min.   :17.50   Min.   : 0.125   no :171  
 1st Qu.: 0.000   male  :286   1st Qu.:27.00   1st Qu.: 4.000   yes:430  
 Median : 0.000                Median :32.00   Median : 7.000            
 Mean   : 1.456                Mean   :32.49   Mean   : 8.178            
 3rd Qu.: 0.000                3rd Qu.:37.00   3rd Qu.:15.000            
 Max.   :12.000                Max.   :57.00   Max.   :15.000            
 religiousness     education       occupation        rating     
 Min.   :1.000   Min.   : 9.00   Min.   :1.000   Min.   :1.000  
 1st Qu.:2.000   1st Qu.:14.00   1st Qu.:3.000   1st Qu.:3.000  
 Median :3.000   Median :16.00   Median :5.000   Median :4.000  
 Mean   :3.116   Mean   :16.17   Mean   :4.195   Mean   :3.932  
 3rd Qu.:4.000   3rd Qu.:18.00   3rd Qu.:6.000   3rd Qu.:5.000  
 Max.   :5.000   Max.   :20.00   Max.   :7.000   Max.   :5.000  
Affairs$utro <- ifelse(Affairs$affairs == 0, 0, 1)#binær variabel

Logistisk regresjon når \(Y\)= 0/1

Vi ønsker å bruke rating og children som prediktorer for utroskap Modellen er \[ log(\frac{p}{1-p})=\beta_0+ \beta_1 rating + \beta_2 children + \epsilon\]

Her er \(p\) sannsynlighet for utroskap, og \(\frac{p}{1-p}\) er oddsen Oddsen er alltid positiv men når vi tar logaritmen så får vi dekket både negative og positive tall!

Tolkningen av \(\beta\)-ene blir mer komplisert!

I R bruker vi glm() funksjonen

log_mod <- glm(utro ~ factor(children) +rating, data=Affairs, family="binomial")
summary(log_mod)

Call:
glm(formula = utro ~ factor(children) + rating, family = "binomial", 
    data = Affairs)

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)          0.28200    0.40264   0.700   0.4837    
factor(children)yes  0.55309    0.24291   2.277   0.0228 *  
rating              -0.47532    0.08608  -5.522 3.36e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 675.38  on 600  degrees of freedom
Residual deviance: 632.80  on 598  degrees of freedom
AIC: 638.8

Number of Fisher Scoring iterations: 4

Tolkning av \(\beta\) i logistisk regresjon

Så barn øker log odds for utroskap, og høyer rating reduserer log odds for utroskap. Vi kan fjerne log ved å bruke eksponentialfunksjon

exp(coef(log_mod))
        (Intercept) factor(children)yes              rating 
          1.3257757           1.7386247           0.6216858

Tolkning: barn øker oddsen for utroskap med 74% og at en enhets økning i rating reduserer oddsen med 38%

hjem