7.5 Dataøving 5

7.5.1 Oppgave 1: Interaktiv øvelse

Før vi tar fatt på dataanalysens begynner vi som vanlig med litt R-trening i swirl. Har du allerede installert pakken swirl (skriv install.packages("swirl") i konsoll hvis ikke) starter du opp swirl med å skrive følgende i konsollen:

library(swirl)
swirl()

Du vil i starten bli bedt om å skrive inn ditt navn og så følger litt info om hvordan swirl fungerer. Du blir så bedt om å velge kurs. Her skal du først velge alternativet ‘R Programming’. Merk at du kanskje må trykke alternativet ‘No. Let me start something new’ for å komme tilbake til hovedmenyen etter å ha brukt swirl tidligere. Du får så se alle modulene dette kurset inneholder. I denne øvingen skal du prøve deg på modul 9 ‘Functions’. I denne modulen vil du lære litt om funksjoner i R.

Merk at det helt til høyre vil står hvor langt du har kommet i prosent. Står du helt fast med et punkt kan du skrive skip() for å hoppe over dette punktet. Når du har fullført en modul blir du spurt om du vil motta ‘credit’ for å ha fullført modulen. Her kan du svare nei. Ønsker du å avbryte underveis skriver du bye(). Skriver du inn det samme navnet når du eventuelt starter swirl igjen kan du fortsette der du slapp. Husk å avslutt swirl (esc) før du begynner på del to av øvingen. Lykke til!

7.5.2 Oppgave 2 - Maskinlæring: Logistisk regresjon og k nærmeste naboer

I denne oppgaven skal vi se på de samme dataene som ble brukt i forelesningen om logistisk regresjon. Vi har data på 10000 kredittkortkunder og vi ønsker å kunne bygge og trene en best mulig modell til å predikere hvilke kunder som vil misligholde sin gjeld.

Oppgave 2.1: Vi starter med å få tak i dataene. Disse er integrert i pakken ISLR. Last inn pakken og ta en titt på dataene ved bruk av følgende linjer (hvordan du kommenterer er opp til deg):

library(ISLR)       # Pakke som inneholder dataene     
head(Default)       # Viser starten på dataframen
str(Default)        # Viser hvilke typer variabler dataframen inneholder
  • Responsvariabelen er
    • : Dette er en kategorisk variabel. Misligholdt kunden gjelden?
  • Forklaringsvariabler:
    • : Kategorisk
    • : Kontinuerlig, størrelsen på gjelden ($)
    • : Kontinuerlig, kundens årlige inntekt ($)

Oppgave 2.2: Det neste vi gjør er å visuelt undersøke avhengigheten mellom det å misligholde (default) og hvor stor gjeld (balance) kunden har . Vanligvis når vi visuelt skal inspisere sammenhengen mellom to variabler lager vi et spredningsplott. Men når den ene variabelen er kategorisk er det mer informativt å sammenligne to boksplott av den kontinuerlige variabelen for hver av gruppene den kategoriske variabelen representer. Dette kan gjøres på følgende måte:

boxplot(balance ~ default, data = Default,
        ylab = "balance", xlab = "default")

Her er det formelen balance ~ default som gjør at boxplot() funksjonen lager to boksplott av balance; et for gruppen som misligholdt (“Yes”) og et for gruppen som ikke misligholdt (“No”). Reflekter over figuren og gjør deg opp en mening om sammenhengen mellom default og balance.

Oppgave 2.3: Det er lurt å dele inn dataene i et treningssett og et testsett når vi driver med maskinlæring. Treningsettet bruker vi til å tilpasse (trene/lære) modellen, mens testsettet bruker vi til å se hvor godt forskjellige modeller presterer. Dette kan gjøres på flere måter, men vi velger her å bruke pakken dplyr som ble beskrevet i siste del av datalabb 2.

Først legger vi til en unik id til hver kunde. Vi lar id-nummeret være lik radnummeret til kunden og til dette bruker vi funksjonen mutate:

library(dplyr)
my_data <- Default %>%
  mutate(id = row_number())

Siden vi nå skal trekke et utvalg av dataene våre kan det være lurt å sikre at resultatet er reprodusibelt ved å sette

set.seed(123)

foran koden som følger. Du kan gjerne velge et annet tall enn 123, men når en gjør dette i forkant av en tilfeldig trekning i R er trekningen bestemt.

Så trekker vi et treningssett bestående av 70 % av dataene ved bruk av funksjonen sample_frac:

train <- my_data %>%  
  sample_frac(.70)

De resterende kundene bruker vi som testsett ved bruk av funksjonen anti_join:

test <- my_data %>%         # Treningssettet er da de resterende 30 % av dataene
  anti_join(train, by = 'id')

Koden over trekker ut alle kunder som ikke har lik id som i treningssettet som derfor svarer til de resterende 30 % av dataene.

Oppgave 2.4: Vi forklarer i denne oppgaven hvordan en logistisk regresjonsmodell kan estimeres, tolkes og brukes. Du vil måtte lage nye modeller med tilsvarende koder i oppgavene som følger.

Vi lager en modell hvor vi bruker variabelen balance (gjeld) som forklaringsvariabel. Vi bruker da funksjonen glm:

model1 <- glm(default ~ balance, data = train,
              family = "binomial")

Syntaksen til glm-funksjonen er veldig lik den vi bruker i regresjon (lm-funksjonen) bortsett fra at vi må spesifisere argumentet family = "binomial" for at å fortelle R at vi ønsker å gjøre en logistisk regresjon. Merk at vi bruker treningssettet til å estimere (trene) modellen ved å spesifisere argumentet data = train.

Det kan være lurt å se om forklaringsvariabelen balance har en signifikant effekt på default ved å bruke summary funksjonen:

summary(model1)

For å tolke hvilken effekt balance (gjeld) har på default (mislighold) er det lurt å regne ut hva effekt en økning på 1 $ i balance har på oddsen for default (Se forelesning):

exp(coef(model1))
##  (Intercept)      balance 
## 2.205746e-05 1.005567e+00

Vi ser at oddsen for default øker med en faktor 1.0056 (en 0.56 % økning) dersom balance øker med 1 $.

Si at du ønsker å predikere sannsynligheten for hvorvidt to kunder med henholdsvis 1000 $ og 2000 $ i balance vil misligholde sitt lån. Da kan vi bruke predict5 på følgende måte:

to_personer <- data.frame(balance = c(1000, 2000))
pred <- predict(model1, newdata = to_personer, type = "response")
pred
##           1           2 
## 0.005648303 0.593966756

Det første argumentet i funksjonen predict er hvilken modell vi skal bruke i prediksjonen (model1). Det andre argumentet newdata er hvilke kunder vi ønsker å predikere misligholdsannsynligheter for. Vi setter dette argumentet til data.frame’n vi har kalt to_personer hvor hver rad svarer til en kunde med et sett forklaringsvariabler (i dette tilfellet to kunder og derfor to rader). Det er viktig at den inneholder en (eller flere) kolonne(r) med kolonnenavn som svarer til navnet til forklaringsvariabelen(e) vi har brukt i modellen. Argumentet type = "response" gjør at vi får returnert sannsynligheten for mislighold og ikke bare verdien av det lineære leddet i modellen.

Hvis vi ut fra disse sannsynlighetene ønsker en klassifiseringsregel som klassifiserer om kunden vil misliholde eller ikke (“Yes/No”) er det naturlig å tildele kunden “Yes” hvis misligholdsannsynligheten overstiger en hvis grense og “No” hvis ikke. Det kan det tenkes at kredittgiver vil være enten konservativ (sette grensen lavt, si 0.3) eller liberal (sette grensen høyt, si 0.7), men i eksempelet under bruker vi en “nøytral” grense på 0.5:

ifelse(pred > 0.5, "Yes", "No")
##     1     2 
##  "No" "Yes"

Som navnet tilsier, vurderer ifelse funksjonen en logisk test i første argumentet (pred > 0.5 hvor pred er misligholdsannsynlighetene vi predikerte) og hvis testen har verdi TRUE gir den ut det du skriver i det andre argumentet ("Yes"), og det tredje argumentet ("No") ellers. Vi ser at individet med 1000 $ i gjeld blir klassifisert som “No”, mens individet med 2000 $ i gjeld blir klassifisert som “Yes”.

Oppgave 2.5: Lag en ny modell ved navn model2 hvor du bruker den kategoriske variabelen student som forklaringsvariabel.

  • Hvilken effekt har det å være student på for oddsen for mislighold?
  • Prediker sannsynligheten for at en student og en ikke-student misligholder gjelden sin.

Oppgave 2.6: Lag en tredje model ved navn model3 hvor du bruker alle forklaringsvariablene.

  • Hvilken effekt har det å være student på for oddsen for mislighold nå? Sammenlign med forrige oppgave.
  • Undersøk visuelt om det er en sammenheng mellom student og balance med å lage et boxplot over balance for studenter og et for ikke-studenter (hint: se Oppgave 2.2).
  • Prediker sannsynligheten for mislighold for en ikke-student og en student med lik balance og income på henholdsvis 1500 $ og 10000 $. Sammenlign med prediksjonen du gjorde i Oppgav 2.4. Basert på model2 og model3, hvordan skal kredittgiver forholde seg til en student versus en ikke-student dersom a) Ingen informasjon om balance eller income er oppgitt og b) dersom en vet balance og income?

Oppgave 2.7: I denne oppgaven skal vi trene opp en knn (k-nærmeste-naboer) modell til å gjøre en tilsvarende klassifisering som den logistiske regresjonsmodellen gjorde over. Funksjonen train som vi trenger er inneholdt i pakken caret (som må installeres ved hjelp av install.packages("caret")). Vi velger å tilpasse en modell hvor antall naboer “k” velges automatisk med kryssvalidering vi argumentet “`trControl”:

library(caret)

# R-kode dersom vi vil velge k automatisk
set.seed(200)
trControl <- trainControl(method  = "cv", # 5-fold kryssvalidering
                          number  = 5)

# Tilpasser modellen
model4 <- train(default ~ balance + income + student,
                 data = train,
                 method = "knn",
                 trControl  = trControl,
                 metric     = "Accuracy")

Vi kan sjekke hvilken k som ble valgt på følgende måte (siden kryssvalidering bruker tilfeldige trekninger kan resultatet bli noe foreskjellig fra gang til gang, selv om datasettet er det samme):

# Hvilken k valgte kryssvalideringen?
k <- model4$finalModel$k
k
## [1] 5

Som for de andre modellene bruker vi funksjonen predict når vi skal predikere og syntaksen er helt lik:

to_kunder <- data.frame(balance = c(1000, 2000), income = 10000, student = c("Yes", "Yes"))
predict(model4, newdata = to_kunder)
## [1] No  Yes
## Levels: No Yes

Merk at i motsetning til de logistiske regresjonsmodellene som predikerte sannsynligheter klassifiserer knn modellen kundene direkte som “Yes”/“No”.

Oppgave 2.8: Vi ønsker å vurdere hvilken av model3 (logistisk regresjon) og model4 (knn) som er best. Vi kan da sjekke hvor godt de klarer å klassifisere testsettet vårt hvor vi vet hvem som har misligholdt lånene sine. Vi starter med å hente ut de sanne verdiene av default i treningssettet:

sann <- test$default  # Den sanne verdien av default i testdataene

Disse skal vi så sammenligne med hvordan modellene klassifiserer de samme kundene basert på de andre variablene. Vi gjør først klassifiseringen med den logistiske regresjonsmodellen:

pred_logreg <- predict(model3, newdata = test, type = "response") # Predikert sannsynlighet
klass_logreg <- ifelse(pred_logreg > 0.5, "Yes", "No")            # Klassifisering av kundene

En oversikt over hvor mange riktige/feil klassifiseringer modellen gjør kan lett oppsummeres med en kontigenstabell. For å lage en kontigenstabell i R bruker vi funksjonen table:

logreg_tab <- table(sann, klass_logreg) # Kontigenstabell
logreg_tab
##      klass_logreg
## sann    No  Yes
##   No  2889   11
##   Yes   70   30

Her kan vi f.eks se at 2889 kunder blir riktig klassifisert som “No”, dvs at de ikke misligholdt lånet og modellen spår at de ikke vil misligholde lånet. Merk at diagonalen (2889 og 30) representerer korrekte klassifiseringer, mens av-diagonal (11 og 70) representer feil klassifiseringer.

Det kan være en fordel å dele kontigenstabellen over med totalt antall kunder for å få andelel i stedet. Funksjonen prop.table gjør nettopp dette:

logreg_tab_norm <- logreg_tab %>% 
  prop.table %>%  # normaliser
  round(3)        # rund av til 3 desimaler
logreg_tab_norm
##      klass_logreg
## sann     No   Yes
##   No  0.963 0.004
##   Yes 0.023 0.010

Summen av diagonalen på denne tabellen (0.963 og 0.01) gir da totalt andel korrekt klassifiseringer:

logreg_tot <- sum(diag(logreg_tab_norm))               # Total andel korrekt klassfisering
logreg_tot
## [1] 0.973

Oppgave 2.9: Gjør en tilsvarende klassifisering av kundene i testsettet med knn modellen model4 og sammenlign med resultatet over. Hvilken modell foretrekker du? Hvilke egenskaper ved klassifisering tror du kredittgiver vektlegger?

7.5.3 Oppgave 3 – Paneldata: Iskaffe og varme dager

I denne oppgaven skal vi jobbe med et simulert paneldatasett. En kaffebarkjede har data fra 150 kaffebarer, observert månedlig over tre år (januar 2021 – desember 2023). For hver kaffebar og måned har vi informasjon om antall iskaffe solgt (iskaffe_salg) og antall dager den måneden med temperatur over 20°C (varme_dager). Kjeden ønsker å vite: Hvor mange ekstra iskaffe selger en typisk bar per ekstra varm dag?

Datasettet er simulert med en kjent, sann effekt – men denne røper vi ikke ennå. Last ned datasettet:

kaffebar_panel.csv

Oppgave 3.1

Les inn datasettet og gjør deg kjent med strukturen:

library(readr)
library(dplyr)

kaffebar <- read_csv("kaffebar_panel.csv")

head(kaffebar)
dim(kaffebar)
  • Hvor mange kaffebarer er det i datasettet?
  • Hvor mange måneder er det observasjoner for?
  • Er panelet balansert (like mange observasjoner per bar)? Bruk gjerne count() og summarise() til å sjekke dette.

Oppgave 3.2

Vi skal nå utforske datasettet visuelt. Begynn med å beregne gjennomsnittlig iskaffe-salg per måned (over alle barer), og lag et linjeplott:

library(ggplot2)

# beregne gjennomsnitt per måned
gj_snitt <- kaffebar %>%
  group_by(maaned_nr) %>%
  summarise(snitt_salg = mean(iskaffe_salg))

# plot
ggplot(gj_snitt, aes(x = maaned_nr, y = snitt_salg)) +
  geom_line() + geom_point() +
  scale_x_continuous(
    breaks = c(1, 7, 13, 19, 25, 31),
    labels = c("Jan 2021", "Juli 2021", "Jan 2022", "Juli 2022", "Jan 2023", "Juli 2023")
  ) +
  labs(x = "Maaned", y = "Gjennomsnittlig iskaffe-salg")

Beskriv hva du ser. Er mønsteret som forventet?

Lag deretter et plott der du velger ut 5 tilfeldige kaffebarer og fargelegger etter kaffebar_id:

set.seed(42)
utvalgte <- sample(unique(kaffebar$kaffebar_id), 5)
kaffebar_utvalg <- kaffebar %>% filter(kaffebar_id %in% utvalgte)

ggplot(kaffebar_utvalg, aes(x = maaned_nr, y = iskaffe_salg, color = factor(kaffebar_id))) +
  geom_line() + geom_point(size = 1) +
  scale_x_continuous(
    breaks = c(1, 7, 13, 19, 25, 31),
    labels = c("Jan 2021", "Juli 2021", "Jan 2022", "Juli 2022", "Jan 2023", "Juli 2023")
  ) +
  labs(x = "Maaned", y = "Iskaffe solgt", color = "Kaffebar")

Hva legger du merke til? Følger barene det samme mønsteret? Ligger de på samme nivå?

Oppgave 3.3

Estimer en enkel regresjonsmodell (pooled OLS, dvs. uten faste effekter) der iskaffe_salg er utfallsvariabel og varme_dager er forklaringsvariabel. Vi bruker feols fra fixest-pakken fordi vi skal bygge videre med faste effekter i neste deloppgave:

library(fixest)

pooled <- feols(iskaffe_salg ~ varme_dager, data = kaffebar)
summary(pooled)
  • Tolk koeffisienten på varme_dager.
  • Tenk tilbake på plottene du laget i 3.2. Kan du se noen grunn til at denne koeffisienten over- eller underestimerer den faktiske effekten av varme dager?

Oppgave 3.4

Estimer nå modellen med faste effekter for kaffebar:

fe <- feols(iskaffe_salg ~ varme_dager | kaffebar_id, data = kaffebar)
etable(pooled, fe)
  • Hva har skjedd med koeffisienten?
  • Hva kontrollerer kaffebar-faste effekter for? Tenk på hva du så i plottet med de 5 utvalgte barene.
  • Hvorfor endrer estimatet seg?

Oppgave 3.5

Legg til tidsfaste effekter (for maaned_nr):

fe_tw <- feols(iskaffe_salg ~ varme_dager | kaffebar_id + maaned_nr, data = kaffebar)
etable(pooled, fe, fe_tw)
  • Sammenlign koeffisientene fra de tre modellene.
  • Hva kan tidsfaste effekter fange opp i denne sammenhengen? Tenk på sesongmønsteret du så i oppgave 3.2.

Oppgave 3.6

Datasettet er simulert med en sann effekt lik \(\beta_1 = 5\): Hver ekstra varm dag gir i gjennomsnitt 5 ekstra iskaffe solgt.

Sammenlign denne sanne verdien med estimatene fra de tre modellene. Hvilken modell kom nærmest, og hvorfor bommer de andre?

7.5.4 Oppgave 4 – Kausal identifikasjon: Ølskatt og trafikkdødsfall

I denne oppgaven bruker vi et ekte datasett med årlige observasjoner for 48 amerikanske delstater fra 1982 til 1988. Datasettet inneholder blant annet informasjon om ølskatt og trafikkdødsfall, og vi ønsker å undersøke sammenhengen mellom disse.

Oppgave 4.1

Installer (om nødvendig) og last inn pakken AER, og last inn datasettet Fatalities:

library(AER)
library(dplyr)
data("Fatalities")

Utforsk datasettet med head(), dim() og str().

  • Hvor mange observasjoner er det?
  • Identifiser variablene state og year som definerer panelstrukturen. Hvor mange delstater og år er det?
  • Er panelet balansert?

Vi skal jobbe med dødsraten per innbygger, ikke det absolutte antallet dødsfall. Variablene fatal (antall trafikkdødsfall) og pop (befolkning) finnes i datasettet. Lag en ny variabel frate som måler antall trafikkdødsfall per 10 000 innbyggere:

Fatalities <- Fatalities %>%
  mutate(frate = fatal / pop * 10000)

Oppgave 4.2

Lag et spredningsplott med beertax (ølskatt i dollar per kasse) på x-aksen og frate (dødsrate) på y-aksen:

library(ggplot2)

ggplot(Fatalities, aes(x = beertax, y = frate)) +
  geom_point() +
  labs(x = "Olskatt (dollar per kasse)", y = "Trafikkdodsfall per 10 000 innb.")

Hva slags sammenheng ser du? Stemmer dette med hva du ville forventet?

Oppgave 4.3

Estimer en pooled OLS-regresjon:

library(fixest)

pooled <- feols(frate ~ beertax, data = Fatalities)
summary(pooled)
  • Tolk koeffisienten. Hva sier den om sammenhengen mellom ølskatt og trafikkdødsfall?
  • Gir fortegnet mening? Hva kan forklare resultatet?

Oppgave 4.4

Estimer modellen med faste effekter for delstat:

fe_state <- feols(frate ~ beertax | state, data = Fatalities)
etable(pooled, fe_state)
  • Hva har skjedd med fortegnet på koeffisienten?
  • Forklar hva delstatsfaste effekter kontrollerer for. Hva slags forskjeller mellom stater kan ha drevet det positive estimatet i pooled OLS?
  • Hva er forskjellen på å sammenligne stater med hverandre (pooled OLS) og å se på endringer innad i samme stat over tid (faste effekter)?

Oppgave 4.5

Legg til tidsfaste effekter:

fe_tw <- feols(frate ~ beertax | state + year, data = Fatalities)
etable(pooled, fe_state, fe_tw)
  • Sammenlign de tre modellene.
  • Hva kan tidsfaste effekter fange opp her? Tenk på nasjonale trender i perioden 1982–1988.

Oppgave 4.6

Refleksjonsoppgave: Er det rimelig å tolke koeffisienten fra modellen med toveis faste effekter som den kausale effekten av ølskatt på trafikkdødsfall? Hva kan eventuelt gjenstå i restleddet som truer en slik tolkning?

  1. Funksjonen predict er satt opp med litt forskjellige argumenter alt ettersom hvilken type modell vi bruker. Du kan lese dokumentasjonen ?predict.glm for å se hvordan den er satt opp for glm objekter↩︎