Modell med én uavhengig variabel

Modeller følger alltid samme oppsett: du velger hvilken modell-funksjon du ønsker å bruke. Inni den presiserer du likningen du skal estimere og dataobjektet hvor du henter informasjonen fra.

Her modellerer jeg respondentens (økonomiske) innvandringsskepsis som en funksjon av vedkommendes arbeidssituasjon.

mod <- lm(formula = Skepsis ~
            Prekaritet,
          data = df)

Som vanlig, lagrer jeg modellresultatet i et eget objekt. Mitt favorittnavn er mod.

Sjekk resultatet ved å be om et sammendrag. Paradoksalt nok, gir summary() mer info enn selve objektet.

#Selve modellobjektet
mod

## 
## Call:
## lm(formula = Skepsis ~ Prekaritet, data = df)
## 
## Coefficients:
## (Intercept)   Prekaritet  
##       4.377        0.430

#Nå sammendraget
summary(mod)

## 
## Call:
## lm(formula = Skepsis ~ Prekaritet, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8066 -1.0432 -0.0432  0.9568  4.9568 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.37657    0.04786  91.450  < 2e-16 ***
## Prekaritet   0.43001    0.10732   4.007 6.53e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.497 on 1220 degrees of freedom
##   (214 observations deleted due to missingness)
## Multiple R-squared:  0.01299,    Adjusted R-squared:  0.01218 
## F-statistic: 16.05 on 1 and 1220 DF,  p-value: 6.527e-05

Modellobjektet er en liste med ymse informasjon knyttet til det.

names(mod)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "na.action"     "xlevels"       "call"          "terms"        
## [13] "model"

Koeffisientene (\(\alpha\) og \(\beta\)) Jeg kan trekke ut koeffisientene fra modellobjektet: De er i praksis en vektor med navn.

koef <- mod$coefficients
koef

## (Intercept)  Prekaritet 
##   4.3765747   0.4300096

Hva er koeffisienten til “Prekaritet” (\(\beta\))?

koef["Prekaritet"]

## Prekaritet 
##  0.4300096

Tekstuell tolkning.

paste("Når prekaritet øker med en skalaenhet, øker innvandringsskepsis med",
      round(koef["Prekaritet"],1),
      "skalaenheter.")

## [1] "Når prekaritet øker med en skalaenhet, øker innvandringsskepsis med 0.4 skalaenheter."

Når prekaritet øker med en skalaenhet, øker innvandringsskepsis med 0.4 skalaenheter. Sagt annerledes: Selv med lik inntekt, vil respondenter med en prekær tilknytning til arbeidslivet være signifikant mer skeptiske til innvandring, selv om effekten er moderat.

Residualene (\(\epsilon_i\)) finner du også her.

mod$residuals[1:10]

##          1          2          3          4          5          6          7 
## -1.3765747  0.2900919  0.9567586  0.6234253 -0.3765747  1.9567586  0.9567586 
##          9         10         11 
##  0.2900919 -0.1399177  2.8600823

Du finner til og med dataene som ble brukt. Her plotter jeg de 6 første observasjonene.

mod$model %>% head()

Modellsammendraget gir et godt overblikk. summary()-funksjonen returner også en liste med statistikk trukket fra modellobjektet.

summary(mod) %>% names

##  [1] "call"          "terms"         "residuals"     "coefficients" 
##  [5] "aliased"       "sigma"         "df"            "r.squared"    
##  [9] "adj.r.squared" "fstatistic"    "cov.unscaled"  "na.action"

Her er spredningen (\(\sigma\)):

summary(mod)$sigma

## [1] 1.497409

Nå kan vi fylle ut likningen.

\[\begin{equation} y_i = \alpha + \beta \times x_i + \sigma^2 \\ y_i = 4.38+ 0.43 \times x_i + 1.5^2 \end{equation}\]

save(mod,
     file = "mod.rda")

Modell med flere uavhengige variabler

Å utvide modellen er uproblematisk. Vanligvis lager jeg da nye modellobjekter, slik at jeg kan ha flere modellresultater for hånd.

Her har jeg kalt den neste modellen min mod2. Her kontrollerer jeg for om respondenten eller hans/hennes foreldre er født utenfor Norge.

mod2 <- lm(formula = Skepsis ~
             Prekaritet +
             Innvandrer,
          data = df)

summary(mod2)

## 
## Call:
## lm(formula = Skepsis ~ Prekaritet + Innvandrer, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.574 -1.085 -0.085  0.915  5.199 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.419      0.051   86.59  < 2e-16 ***
## Prekaritet     0.439      0.107    4.10  4.5e-05 ***
## Innvandrer    -0.284      0.118   -2.41    0.016 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.5 on 1217 degrees of freedom
##   (216 observations deleted due to missingness)
## Multiple R-squared:  0.0177, Adjusted R-squared:  0.0161 
## F-statistic:   11 on 2 and 1217 DF,  p-value: 1.89e-05

Modellsammendraget gjør det klart at respondenter med innvandrerbakgrunn er mindre innvandringsskeptiske.

Din tur!

• Kan du lage en ny modell hvor du kontrollerer for respondentens objektive inntekt (“Inntekt”) og subjektive oppfattelse av egen inntekt (“SubjektivInntekt”)? Kall det nye modellobjektet mod3.
• \(R^2\) er rapportert i modellsammendraget. Hvor mye øker \(R^2\) med når du legger til disse variablene?
• sammenstill modellene ved hjelp av stargazer(mod, mod2, ...).
• lek deg rundt med tilleggsargumenter i stargazer().

Her er mitt forsøk.

#Ny modell
mod3 <- lm(Skepsis ~ 
             Prekaritet + 
             Innvandrer +
             Inntekt,
           data = df)
#Resulattabell
stargazer(mod, mod2, mod3,
          type = "html",
          #Tittel
          title = "Arbeidslivstilknytningens pavirkning pa innvandringsskepsis",
          #Vektor med variabelnavn + konstantledd
          covariate.labels = c("Usikker tilknytning til arbeidslivet",
                               "Innvandrer",
                               "Inntekt",
                               "Konstantledd"),
          #Hva er min avhengige variabel?
          dep.var.labels = "Innvandringsskepis",
          #Endre til norsk kolonnenavn
          dep.var.caption = "Avhengig variabel:",
          #Komma som desimalskilletegn
          decimal.mark = ",",
          #Rapportere hva slags modell er det i objektet?
          model.names = TRUE)

Koeffisientplot med en kategorisk, ikke-kontinuerlig x-variabel

nyedata <- 
  #lag datasett
  data.frame(
    #presiser variabelnavn og verdier; korte variabler "resirkuleres"
  Prekaritet = c(0, 1))

#ta en titt
nyedata

predict(mod, 
        newdata = nyedata)

##   1   2 
## 4.4 4.8

predict(mod, 
        newdata = nyedata,
        se = TRUE)

## $fit
##   1   2 
## 4.4 4.8 
## 
## $se.fit
##     1     2 
## 0.048 0.096 
## 
## $df
## [1] 1220
## 
## $residual.scale
## [1] 1.5

predikerte <- predict(mod, 
                      newdata = nyedata,
                      se = TRUE)

#Opprett et nytt datasett
dfp = nyedata

#legg til de predikerte verdiene
dfp$estimat <- predikerte$fit 
dfp$hoymargin <- predikerte$fit + predikerte$se.fit
dfp$lavmargin <- predikerte$fit - predikerte$se.fit

dfp

#hent pakken fra biblioteket
library(ggplot2)

#opprett et tomt ark
ggplot(dfp) +
  #legg til punkter
  geom_point(aes(
    #De predikerte y-verdiene
    y = estimat,
    #Ulike scenarier
    x = Prekaritet))

library(ggplot2)

ggplot(dfp) +
  #Plott punktene
  geom_point(aes(
    #De predikerte y-verdiene
    y = estimat,
    #Ulike scenarier
    x = Prekaritet)
  ) +
  #Feilbarrer
  geom_errorbar(aes(
    #Scenariene
    x = Prekaritet,
    #Hvor søylen begynner
    ymin = lavmargin,
    #Hvor søylen ender
    ymax = hoymargin)
  )

Scenarier på speed: kontinuerlige variabler

Når \(x\)-variaben vår er kontinuerlig, kan vi illustrere effekten som en regresjonslinje.

Dette er samme logikken som med to scenarier, men denne gangen lar vi inntekt gå fra lav til høy på kontinuerlig vis.

Jeg estimerer modellen min.

mod1 <- lm(Skepsis ~
             Inntekt +
             Prekaritet,
           data = df)

seq(from = 0, 
    to = 10,
    by = 0.1)

nyedata <- data.frame(
  #hold prekaritet konstant
  Prekaritet = c(0),
  #la inntekt gå fra 0 til 10 
  Inntekt = seq(from = 0,
                to = 10, 
                by = 0.1))
head(nyedata)

predikerte <- predict(mod1, 
                      newdata = nyedata,
                      se = TRUE)
#Opprett et nytt datasett
dfp = nyedata

dfp$estimat <- predikerte$fit 
dfp$hoymargin <- predikerte$fit + predikerte$se.fit
dfp$lavmargin <- predikerte$fit - predikerte$se.fit

ggplot(dfp) +
  geom_line(aes(
    y = estimat,
    x = Inntekt)
  )

ggplot(dfp) +
  geom_line(aes(
    y = estimat,
    x = Inntekt)
  ) +
  geom_ribbon(aes(ymax = hoymargin,
                  ymin = lavmargin,
                  x = Inntekt),
              #gjennomsiktig
              alpha = 0.5)

ggplot(dfp) +
  geom_line(aes(
    y = estimat,
    x = Inntekt)
  ) +
  #hoyeste grenseverdi
  geom_line(aes(
    y = hoymargin,
    x = Inntekt),
    lty = 2 #stipulert linje
  ) +
  #laveste grenseverdi
  geom_line(aes(
    y = lavmargin,
    x = Inntekt),
    lty = 2 #stipulert linje
  )

Scenarier for de viderekommende: Simulere utfallene

Fram til nå har vi antatt at feilmarginene er de samme uansett hvor i prediksjonen du er og at det beste estimatet vår er en gjennomsnittsverdi. Det er ikke alltid tilfelle; særlig når du går over til GLM-er (f.eks. logit). Da kan vi simulere utfall.

Jeg gjør simuleringen i MASS-pakken, så jeg begynner med å hente den inn fra biblioteket. Om du ikke har brukt den før, må du innstallere den først (install.packages("MASS"))

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

Pakken inneholder uheldigvis enkelte funksjoner som har samme navn som funksjoner i dplyr. Dette får du beskjed om når du laster pakken inn. Her er det select() som blir overskrevet av MASS. Om du likevel ønsker å bruke select-funksjonen fra dplyr, kan du presisere dette i koden din dplyr::select().

#Sjekk variablene i modellen
names(mod$coefficients)

## [1] "(Intercept)" "Prekaritet"

#Lag scenarier: her en prekær og en ikke-prekær arbeidstaker
nyedata <- data.frame("Intercept" = 1, #Her må jeg presisere konstantleddet også
                      "Prekaritet" = c(0,1))

#Simuler
sim <-MASS::mvrnorm(
  #Antall simuleringer
  n = 10000,
  #Koeffisientene
  mu= mod$coefficients,
  #Usikkerheten/spredningen
  Sigma= vcov(mod)
)

#Multipliser nye koeffisienter med de hypotetiske dataene
prediksjoner <- as.matrix(nyedata) %*% t(sim)

#Hver linje i datasettet er nå et scenario.
dim(prediksjoner)

## [1]     2 10000

#... og en medianverdi
median(prediksjoner[1,])

## [1] 4.4

median(prediksjoner[2,])

## [1] 4.8

# bruk apply() for å beregne kvantilene fra simuleringen din
kvantiler <- apply(X = prediksjoner,
                   #Margin (i == linje)
                   MARGIN = 1,
                   #Funksjonen som skal brukes på denne linjen
                   FUN = quantile,
                   #Tilleggsargument
                   na.rm = T,
                   #Enda et: 95% konfidensintervall
                   probs = c(.025,.5,.975))

#Tilt dataene slik at linje == scenario
dfp <- as.data.frame(t(kvantiler))

Legg de hypotetiske dataene til dfp.

dfp <- cbind(dfp,
             nyedata)

dfp

Plott i vei!

library(ggplot2)

ggplot(dfp) +
  geom_point(aes(
    #De predikerte y-verdiene
    y = `50%`,
    #Ulike scenarier
    x = Prekaritet))

library(ggplot2)

dfp$Arbeidslivstilknytning <- c("Trygg",
                                "Prekaer")

ggplot(dfp) +
  #Plott punktene
  geom_point(aes(
    #De predikerte y-verdiene
    y = `50%`,
    #Ulike scenarier
    x = Prekaritet,
    #Legg til en kategorisk variabel
    color = Arbeidslivstilknytning)
  ) +
  #Feilbarrer
  geom_errorbar(aes(
    #Scenariene
    x = Prekaritet,
    #Hvor søylen begynner
    ymin = `2.5%`,
    #Hvor søylen ender
    ymax = `97.5%`,
    #Fjern hakene
    width = 0,
    #grupper etter kategorisk variabel
    color = Arbeidslivstilknytning)
  ) +
  xlim(-.25, 1.25) +
  #Skaler
  ylim(3,6) +
  scale_color_manual(values = c("grey", "black")) +
  #bruker et annet "tema" som basis for estetikken
  theme_light() +
  #... men endrer innstillingene
  #Siden x er kategoris, fjerner jeg x-verdiene
  theme(axis.text.x= element_blank(),
        axis.ticks.x = element_blank(),
        #Fjern vertikale hjelpelinjer
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        #Hvor skal tegnforklaringen vaere
        legend.position = "bottom",
        #Retning
        legend.direction = "horizontal",
        legend.title = element_text(face = "bold",
                                    size = 10)) +
  ggtitle(label = "Effekten av arbeidslivstilknytning pa innvandringsskepsis",
          subtitle = "ESS (2014)") +
  ylab("Innvandringsskepsis") +
  xlab("")

Regn ut de simulerte grenseverdiene: det mest sannsynlige utfallet (medianen) og grenseverdiene for prediksjonsintervallet (95% mest sannsynlige utfall).

dfp <-
prediksjoner %>%
  #Transposer datasettet slik at hvert scenario er en variabel
  t() %>%
  #Sørg for at dette er et datasett
  as.data.frame %>%
  #Lag sammendrag
  summarize(
    #Funksjonen across åpner for å gjøre den samme operasjonen på mange datasett
    across(
      #Variablene som skal sammendras
      1:2,
      #Kvantilene (95% intervall)
           quantile,
                   c(0.025, 0.5, 0.975))) %>%
  #Snu datasettet igjen
  t() %>%
  as.data.frame() %>%
  #Legg til info om sceariene
  mutate(x = 0:1,
         Arbeidslivstilknytning = c("Trygg", "Prekær"))

names(dfp)[1:3] <- c("ymin", "y", "ymax")

library(ggplot2)

ggplot(dfp) +
  #Plott punktene
  geom_point(aes(
    #De predikerte y-verdiene
    y = y,
    #Ulike scenarier
    x = x)
  ) +
  #Feilbarrer
  geom_errorbar(aes(
    #Scenariene
    x = x,
    #Hvor søylen begynner
    ymin = ymin,
    #Hvor søylen ender
    ymax = ymax,
    #Fjern hakene
    width = 0)
  )  +
  #Skaler
  # ylim(4,5) +
  coord_cartesian(xlim = c(-0.5, 1.5),
                  ylim = c(4,5))+
  scale_x_continuous(breaks=c(0, 1),
                     labels=c("Fast arbeid",
                              "Midlertidig/uten \narbeid")) +
  #bruker et annet "tema" som basis for estetikken
  theme_minimal() +
  #... men endrer innstillingene
  #Siden x er kategoris, fjerner jeg x-verdiene
  theme(#axis.text.x= element_text(angle = 45),
        axis.ticks.x = element_blank(),
        axis.text.x = element_text(face = "italic"),
        #Fjern vertikale hjelpelinjer
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()) +
  ggtitle(label = "Effekten av arbeidslivstilknytning pa innvandringsskepsis",
          subtitle = "ESS (2014)") +
  ylab("Innvandringsskepsis (0-10)") +
  xlab("Arbeidslivstilknytning")