MET4 HEKS løsningsforslag V24

Oppgave 1: Deskriptiv statistikk og utforskende dataanalyse

a) Utfør en utforskende dataanalyse med deskriptiv statistikk basert på det balanserte datasettet whr_balansert. Fokuser spesielt på variabelen lifeladder. Presenter en tabell og tre figurer dere mener er relevante, og forklar hva de viser.

Det er flere ulike måter å utforske et datasett på så fasit for denne oppgaven er bare et eksempel. Noe det likevel kan være lurt å fokusere på i en slik utforskende analyse er om det er manglende verdier for noen variabler (NAs), distribusjon av variablen man er interessert i, og om det er betydelig variasjon på tvers av grupper.

# Tabell 1

sumstats <- whr_balansert %>%
  select(lifeladder, logGDP, socialsupport, lifeexpectancy, 
         freedom, percepcorruption, generosity)

library(stargazer)
stargazer(sumstats, type="text") 

==================================================
Statistic         N   Mean  St. Dev.  Min    Max  
--------------------------------------------------
lifeladder       960 5.782   1.023   2.694  7.889 
logGDP           868 9.780   0.924   7.586  11.371
socialsupport    879 0.841   0.103   0.434  0.987 
lifeexpectancy   860 66.516  5.649   48.120 75.200
freedom          873 0.775   0.125   0.333  0.964 
percepcorruption 832 0.743   0.189   0.151  0.977 
generosity       861 -0.016  0.156   -0.335 0.553 
--------------------------------------------------

Fra denne tabellen ser vi at snittverdien for lifeladder er 5.782, og det er 960 observasjoner i denne variabelen (N = 960). På de andre variablene derimot ser vi at det er en del manglende observasjoner (N < 960, kodet som NA). Det betyr at når vi senere skal gjennomføre regresjonsanalyser, vil vi ikke ha informasjon om kontrollvariablene for alle observasjonene der vi har informasjon om utfallsvariabelen lifeladder.

# Histogram for Life Ladder
ggplot(whr_balansert, aes(x = lifeladder)) +
  geom_histogram(binwidth = 0.1) +
  theme_minimal() +
  labs(title = "Histogram of Life Ladder Scores", x = "Life Ladder Score", y = "Frequency")

# Boksplott for Life Ladder fordelt på region
# Kan også bruke boxplot(lifeladder ~ region, data = whr_balansert)
# I ggplot:
ggplot(whr_balansert, aes(x = region, y = lifeladder)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Life Ladder Scores by Region", x = "Region", y = "Life Ladder Score") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Scatterplot for sammenhengen mellom lykke og Log BNP
ggplot(whr_balansert, aes(x = logGDP, y = lifeladder)) +
  geom_point() +
  geom_smooth(method='lm', formula= y~x)

Historgramet viser hvor ofte ulike lifeladder nivår blir rapportert, og vi ser at det er en nokså tilnærmet normalfordeling med en liten opphopning rundt 7. Fra boxplottet ser vi at det er noe regional variasjon i lykkenivåer, men det er kanskje overraskende stabilt hvis vi ser vekk fra Nord Amerika og Australia/New Zealand og Sør-Asia. Likevel ser vi stor variasjon i maks- og minimumsverdiene i de fleste regioner, som indikerer at det er noen land-år observasjoner som avviker fra snittet i regionen. Videre ser vi fra scatterplottet at det er en linær og positiv korrelasjon mellom lykkenivå og log BNP.

b) Koronaåret 2020 var et år mange ser tilbake på som tøft. Finner vi grunnlag for dette i datasettet whr_balansert? Formuler en hypotese og test den.

En måte å teste om 2020 var et ulykkelig år er å anta at folk var like lykkelige i 2020 og 2019, og teste den hypotesen mot alternativhypotesen om at folk var mindre lykkelige i 2020. Til dette bruker vi en t-test for to populasjoner og en ensidig test pga. antagelsen vår om at 2020 var et tøft år.

Vi ser fra denne hypotesetesten at 2020 ikke ser ut til å være statistisk sett mindre lykkelig enn 2019, med en p-verdi på 0.4921. Derfor kan vi ikke forkaste null-hypotesen om at folk var like lykkelig i 2020 som i 2019.

# Var 2020 et ulykkelig år? 
whr_2020 <- whr_balansert %>%
  filter(year == 2020)
whr_2019 <- whr_balansert %>%
  filter(year==2019)

# Test om lykken i 2020 var lavere enn lykken i 2019, med en ensidig t-test.
ttest_2020 <- t.test(whr_2020$lifeladder, whr_2019$lifeladder, alternative = "less")
print(ttest_2020)

    Welch Two Sample t-test

data:  whr_2020$lifeladder and whr_2019$lifeladder
t = -0.019873, df = 156.84, p-value = 0.4921
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
      -Inf 0.2580856
sample estimates:
mean of x mean of y 
 5.863012  5.866150 

Alternativt kan man sette opp en paret t-test som tar høyde for forskjeller mellom land fordi målinger av lykke i samme land i 2019 og i 2020 vil med stor sannsynlighet være avhengige.

Dette kan være en fordel i en situasjon der man enten ikke har perfekt randomisering av behandling (“treatment”), og man er redd det er forskjeller mellom landene man ser på som ikke skyldes at de var utsatt for ulik behandling (“treatment”). Det kan også være en fordel å bruke paret t-test hvis man har mye variasjon i datasettet, der noen land er veldig lykkelige og andre er mindre lykkelige.

I dette tilfellet har vi måling for alle land før korona og med korona i whr_balansert, og vi ønsker dermed kun å ta høyde om det er for stor variasjon mellom landene. Fra testen under ser vi at dette ikke hadde påvirkning på konklusjonen, hvor en paret t-test gir en p-verdi på 0.946 og man kan ikke forkaste null-hypotesen om at 2019 og 2020 var like lykkelige.

# Paret T-test hvor vi bruker matchede par av land og tester differeansen mellom 2020 og 2019 (regn ut for hånd eller bruk R)
D <- whr_2020$lifeladder - whr_2019$lifeladder 
gjennomsnitt_D <- mean(D)
st_avvik_D <- sd(D)

t.test(D, mu = 0)

    One Sample t-test

data:  D
t = -0.067109, df = 79, p-value = 0.9467
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.09619531  0.08992031
sample estimates:
 mean of x 
-0.0031375 
t.test(whr_2020$lifeladder, whr_2019$lifeladder, paired =  TRUE)

    Paired t-test

data:  whr_2020$lifeladder and whr_2019$lifeladder
t = -0.067109, df = 79, p-value = 0.9467
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -0.09619531  0.08992031
sample estimates:
mean difference 
     -0.0031375 

Men husk fra introduksjonen at det manglet mange land fra 2020, som er blitt helt ekskludert fra datasettet whr_balansert. Hvis det å være observert i datasettet korrelerer med lykkenivå i landet, vil vi potensielt trekke feil konklusjon og vi kan derfor ikke per nå generalisere utover det sampelet vi har.

Alt i alt betyr det at det kanskje var lavere (eller høyere) lykkenivå hos de landene vi ikke har data på, og det kan også ha vært noen grupper innad i land med lavere lykke i 2020, men som ikke har påvirket landsgjennomsnittet nok til å gi utslag i denne testen. Dette kan skje, for eksempel, ved at også noen grupper ble mer lykkeligere og dermed endret ikke landsgjennomsnittet seg.

Oppgave 2: Vurdering av populasjon

a) Bruk deskriptiv statistikk basert på datasettet whr_average til å utforske om filtreringen vi gjør for å lage det balanserte datasettet whr_balansert er tilfeldig. Fokuser spesielt på variablene lifeladder og logGDP.

ggplot(whr_average, aes(x = included_in_balanced, y = lifeladder)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Life Ladder Scores by filtration", x = "Included in balanced?", y = "Life Ladder Score") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(whr_average, aes(x = included_in_balanced, y = logGDP)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Log-GDP by filtration", x = "Included in balanced?", y = "Log-GDP") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
Warning: Removed 6 rows containing non-finite values (`stat_boxplot()`).

Fra boxplottene kan det set ut til at land som inkluderes i whr_balansert er lykkeligere og rikere enn landene som eksluderes.

b) Estimer en eller flere logistiske regresjonsmodeller som kan brukes til å utforske om filtreringen vi gjør for å lage det balanserte datasettet whr_balansert er tilfeldig. Diskuter hvilken populasjon whr_balansert eventuelt representerer.

Dersom det er helt tilfeldig hvilke land som inkluderes i whr_balansert så skal \(P(\texttt{included{\_}in{\_}balanced} = \texttt{TRUE})\) være uavhengig av forklaringsvariablene. For å undersøke dette tilpasser vi modellene:

logitmod1 <- glm(included_in_balanced ~ lifeladder + logGDP + socialsupport + lifeexpectancy + freedom + percepcorruption + generosity,
                data = whr_average,
                family = "binomial")

logitmod2 <- glm(included_in_balanced ~ lifeladder + logGDP + socialsupport + lifeexpectancy + freedom + percepcorruption + generosity + region,
                data = whr_average,
                family = "binomial")

stargazer(logitmod1, logitmod2, type = "text")

=====================================================================
                                             Dependent variable:     
                                         ----------------------------
                                             included_in_balanced    
                                              (1)            (2)     
---------------------------------------------------------------------
lifeladder                                   0.166          0.029    
                                            (0.398)        (0.484)   
                                                                     
logGDP                                       0.212         -0.027    
                                            (0.362)        (0.402)   
                                                                     
socialsupport                                0.216         -0.636    
                                            (2.815)        (2.914)   
                                                                     
lifeexpectancy                               0.098*         0.051    
                                            (0.053)        (0.069)   
                                                                     
freedom                                      0.565          3.330    
                                            (2.092)        (2.585)   
                                                                     
percepcorruption                             2.475*        3.649**   
                                            (1.279)        (1.699)   
                                                                     
generosity                                   -1.058        -1.970    
                                            (1.383)        (1.718)   
                                                                     
regionCommonwealth of Independent States                   -1.214    
                                                           (1.024)   
                                                                     
regionEast Asia                                             0.030    
                                                           (1.425)   
                                                                     
regionLatin America and Caribbean                         -2.095**   
                                                           (0.952)   
                                                                     
regionMiddle East and North Africa                         -1.640*   
                                                           (0.945)   
                                                                     
regionNorth America and ANZ                                15.837    
                                                         (1,174.636) 
                                                                     
regionSouth Asia                                           -2.421*   
                                                           (1.281)   
                                                                     
regionSoutheast Asia                                       -2.194*   
                                                           (1.232)   
                                                                     
regionSub-Saharan Africa                                  -2.295**   
                                                           (1.167)   
                                                                     
regionWestern Europe                                       -0.019    
                                                           (1.157)   
                                                                     
Constant                                   -11.549***      -6.408    
                                            (2.847)        (4.535)   
                                                                     
---------------------------------------------------------------------
Observations                                  159            159     
Log Likelihood                              -92.592        -84.646   
Akaike Inf. Crit.                           201.183        203.293   
=====================================================================
Note:                                     *p<0.1; **p<0.05; ***p<0.01

Tilfeldig hvilke land som blir inkludert?: For å unngå problemet med multippel testing burde vi ideelt sett ha testet \(\beta_1=\beta_2=...=0\) som i vanlig lineær regresjon svarer til en F-test. Vi har ikke lært om en tilsvarende test for logistisk regresjon i MET4. Vi ser allikevel at flere (enn det vi forventer) av forklaringsvariablene er signifikante i begge modellene. Fra modell (1) ser vi at for land som inkluderes er det lengre forventet levetid og mer korrupsjon. Selv om de ikke er signifikant forskjellig fra null noterer vi at verdien av koeffisientene til lifeladder og logGDP tilsier at de inkluderte landene er rikere og lykkeligere. Fra modell (2) der vi også inkluderer verdensdel ser vi at land som blir inkludert har signifikant lavere sjanse å komme fra Latin-Amerika, Karibien, Afrika og Asia sammenlignet med land som ikke blir inkludert. Vi ser også at når vi kontrollerer for verdensdel er ikke levetid lenger signifikant, samt at størrelsesorden på koeffsienten til lifeladder og logGDP omtrent er null. Dette skyldes at verdensdel og disse variablene er avhengige: (Indikasjon på at hvis du er fra en av de nevnte verdensdelene har du i gjennomsnitt lavere forventet levetid, er fattigere og mindre lykkeligere).

Hvilken populasjon representer whr_balansert?: Ut fra den logistiske regresjonsmodellen vil datasettet whr_balansert representere en populasjon med større innslag av land hvor det er lang levetid og mer korrupsjon. (Dette med korrupsjon virker med første øyekast litt rart, se oppgave 3a for sammenhengen mellom lykke og korrupsjon). I praksis betyr dette at land fra verdensdelene Latin-Amerika og Karibien, Afrika og Asia vil også være underrepresentert i dette datasettet.

c) Bruk først hele datasettet whr til å se på gjennomsnitt av logGDP og lifeladder over tid. Er det noen overraskende funn? Sammenlign med tilsvarende analyse for det balanserte datasettet whr_balansert. Kommenter i lys av funnene dine i oppgave 2a) og 2b).

Vi begynner med å plotte trenden for lifeladder og logGDP fra datasettet whr:

whr_average_per_year <- whr %>%
  group_by(year) %>%
  summarise(average_lifeLadder = mean(lifeladder, na.rm = TRUE),
            average_logGDP = mean(logGDP, na.rm = TRUE)
            )

ggplot(whr_average_per_year, aes(x = year, y = average_lifeLadder)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Time Trend of Average Life Ladder",
       x = "Year",
       y = "Average Life Ladder") +
  scale_x_continuous(breaks=c(2012, 2014, 2016, 2018, 2020, 2022))

ggplot(whr_average_per_year, aes(x = year, y = average_logGDP)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Time Trend of Average Log GDP per Capita",
       x = "Year",
       y = "Average Log GDP per Capita") +
  scale_x_continuous(breaks=c(2012, 2014, 2016, 2018, 2020, 2022))

Overaskende funn: Folk ser veldig lykkelig ut i koronaåret 2020 - suspekt. Den andre figuren viser høy BNP per capita i 2020, et år vi vet det ikke gikk så bra økonomisk på verdensbasis - også suspekt.

Kommentar i lys av 2a) og 2b) Fra oppgave 2a) og 2b) vet vi at landene som ble filtrert bort i gjennomsnitt er fattigere og mer ulykkeligere (dog ikke signfikant). Vi får oppgitt i introduksjonen at mange land blir filtrert bort fordi vi mangler observasjoner for dem i 2020 i originaldataene whr. Det er derimot rimelig å tro at vi har observasjoner for mange av disse landene i 2019 og 2021 i det samme datasettet. Så når disse landene er inkludert i gjennomsnittene for 2019 og 2021, men ikke i 2020 i figurene over så får vi disse kunstige toppene i 2020.

Vi tar så en titt på de tilsvarende figurene i det balanserte datasettet:

whr_average_per_year <- whr_balansert %>%
  group_by(year) %>%
  summarise(average_lifeLadder = mean(lifeladder, na.rm = TRUE),
            average_logGDP = mean(logGDP, na.rm = TRUE)
            )
# Endre skala her, slik at de er sammenlignbare med de over, på en fast skala. Lettere å sammenligne og vil gjøre at det ikke fremstår som en så ekstrem endring. 
ggplot(whr_average_per_year, aes(x = year, y = average_lifeLadder)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Time Trend of Average Life Ladder",
       x = "Year",
       y = "Average Life Ladder") +
  scale_x_continuous(breaks = c(2012, 2014, 2016, 2018, 2020, 2022)) +
  scale_y_continuous(limits = c(5.4, 6))  # Sett y-aksen liknende tidligere figur for å lettere sammenligne 

ggplot(whr_average_per_year, aes(x = year, y = average_logGDP)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Time Trend of Average Log GDP per Capita",
       x = "Year",
       y = "Average Log GDP per Capita") +
  scale_x_continuous(breaks=c(2012, 2014, 2016, 2018, 2020, 2022)) +
  scale_y_continuous(limits = c(9,10)) # Sett y-aksen liknende tidligere figur for å lettere sammenligne 

Disse figurene er kanskje mer som forventet. Vi ser en dupp i logGDP i 2020, og vi har ikke en lykketopp i 2020, i samsvar med konklusjonen i 1b).

Oppgave 3: Regresjonsanalyse

a) Utfør en eller flere regresjonsanalyser for å undersøke forholdet mellom lifeladder (som avhengig variabel) og ulike forklaringsvariabler. Diskuter resultatenes tolkning i lys av antagelsen om eksogenitet. Bruk datasettet whr_balansert.

Før vi setter igang med regresjonsanalysene bør vi sjekke avhengighetsstrukturen mellom variablene, hvis ikke dette ble gjort i oppgave 1a. Dette kan vi gjøre ved å regne ut korrelasjonen mellom ulike potensielle kontrollvariabler og utfallsvariabelen lifeladder, og deretter plotte dem for å se etter ikke-linearieter. Om det ikke er linære sammenhenger vil det indikere at vi bør teste ut log-transformering av variabelene før vi inkluderer dem i regresjonsanalysen.

whr_balansert %>%
  select(-year, -region, -country) %>% # Tar vekk ikke-numeriske variabler 
  cor(use = "complete.obs")  # Ekskluder NAs 
                 lifeladder      logGDP socialsupport lifeexpectancy    freedom
lifeladder        1.0000000  0.78857847     0.6894389     0.73929073  0.4631776
logGDP            0.7885785  1.00000000     0.6712191     0.87569480  0.1982224
socialsupport     0.6894389  0.67121913     1.0000000     0.55962339  0.2757587
lifeexpectancy    0.7392907  0.87569480     0.5596234     1.00000000  0.2117747
freedom           0.4631776  0.19822245     0.2757587     0.21177472  1.0000000
percepcorruption -0.5800848 -0.43421173    -0.2626353    -0.40317889 -0.5545722
generosity        0.2509389  0.04127431     0.1242181     0.05790208  0.4185880
                 percepcorruption  generosity
lifeladder             -0.5800848  0.25093894
logGDP                 -0.4342117  0.04127431
socialsupport          -0.2626353  0.12421813
lifeexpectancy         -0.4031789  0.05790208
freedom                -0.5545722  0.41858798
percepcorruption        1.0000000 -0.35244414
generosity             -0.3524441  1.00000000

Her ser vi at lifeladder er sterkt positivt korrelert med BNP, sosial støtte, og forventet levealder. Det er også en ganske sterk negativ korrelasjon mellom lykke og oppfattet korrupsjonsnivå.

# Parvise plot 
pairs(lifeladder~., data = whr_balansert %>% select(-c(year, region, country)))

Fra første linje ser vi plott av lifeladder plottet mot de ulike variablene. Her ser vi at det stort sett er linære sammenhenger, men vi kan se fra freedom, percepcorruption, og generosity variabelene at det ser ut til å være noe hetroskedastisitet.

For percepcorruption ser vi noe som tyder på en ikke-linær sammenheng, så vi kan også teste den som log-transformert.

Det er forøvrig en litt artig sammenheng mellom korrupsjon og lykke: De som har lav precorruption har en meget stabil høy lykke, men så snart precorruption bikker rundt 0.6 så er det stor variasjon i lykke. Sagt på en annen måte, i lite korrupte land er folk lykkelige, mens du kan både være lykkelig og ulykkelig i korrupte land!

Kommentar: En klarer aldri å transformere seg vekk fra heteroskedastisitet, da må man modellere variansen som funksjon av precorruption også og det har vi ikke lært i MET4.

pairs(log(lifeladder) ~ log(percepcorruption), data = whr_balansert)

pairs(lifeladder ~ percepcorruption, data = whr_balansert)

Fra disse plottene ser vi at det er uten log-transformering som fremstår som mest linær, og vi fortsetter å benytte variabelen uten log-transformasjon senere i oppgaven. Merk at selv om vi har brukt pairs() i MET4, så er det også helt greit å se på ett og ett spredningsplott.

Da kan vi gjennomføre noen regresjonsanalyser. Først tester vi en enkel regresjons-modell (OLS) uten faste effekter. Men siden vi vet at det er en del variasjon mellom land (som vi så i boksplottet til regioner hvor det var en del variasjon og outliers) og utvikling over tid, så utvider vi modellen vår til å inkludere først tid- og region-faste effekter og til slutt tid- og land-faste effekter.

library(stargazer)
# Lineær regresjonsmodell
lm_model <- lm(lifeladder ~ logGDP + socialsupport + freedom + percepcorruption +
                  lifeexpectancy + generosity, data = whr_balansert)

# Inkluder region- og tidsfaste effekter
fe_model_region <- lm(lifeladder ~ logGDP + socialsupport + freedom + percepcorruption
                      + lifeexpectancy + generosity + region + as.factor(year),
                      data = whr_balansert)


# Inkluder land- og tidsfaste effekter 
fe_model_twfe <- lm(lifeladder ~ logGDP + socialsupport + freedom + percepcorruption 
                    + lifeexpectancy + generosity + country + as.factor(year),
                    data = whr_balansert)

stargazer(lm_model, fe_model_region, fe_model_twfe, type = "text", omit = c("country", "region", "year"))

===============================================================================================
                                                Dependent variable:                            
                    ---------------------------------------------------------------------------
                                                    lifeladder                                 
                              (1)                       (2)                      (3)           
-----------------------------------------------------------------------------------------------
logGDP                      0.380***                 0.515***                  1.748***        
                            (0.046)                   (0.047)                  (0.185)         
                                                                                               
socialsupport               2.546***                 2.115***                  1.998***        
                            (0.234)                   (0.237)                  (0.325)         
                                                                                               
freedom                     1.330***                 0.924***                   0.366*         
                            (0.181)                   (0.197)                  (0.218)         
                                                                                               
percepcorruption           -0.974***                 -0.866***                -1.009***        
                            (0.125)                   (0.128)                  (0.221)         
                                                                                               
lifeexpectancy              0.034***                 0.028***                 -0.101***        
                            (0.006)                   (0.009)                  (0.023)         
                                                                                               
generosity                  0.432***                 0.777***                  0.542***        
                            (0.127)                   (0.137)                  (0.184)         
                                                                                               
Constant                   -2.666***                 -3.058***                 -5.275**        
                            (0.323)                   (0.581)                  (2.193)         
                                                                                               
-----------------------------------------------------------------------------------------------
Observations                  797                       797                      797           
R2                           0.778                     0.814                    0.924          
Adjusted R2                  0.776                     0.808                    0.914          
Residual Std. Error     0.495 (df = 790)         0.459 (df = 771)          0.308 (df = 701)    
F Statistic         460.636*** (df = 6; 790) 134.838*** (df = 25; 771) 89.648*** (df = 95; 701)
===============================================================================================
Note:                                                               *p<0.1; **p<0.05; ***p<0.01

Fra disse 3 linære modellene ser vi at sammenhengen mellom log BNP og lykke er statistisk signifikant forskjellig fra 0 på et 1% signifikansnivå, men punktestimatet varierer en god del mellom modellene. For modell (1) er punktestimatet til b1 lik 0.380 som sier oss at ved en 1%-økning i BNP er det assosiert med 0.380/100=0.0038 enhets økning i lykkenivå. Det at punktestimatet varierer mye mellom modellene forteller oss at det trolig er endogenitet i modellen, og at b1 fra modell (1) ikke er et forventningsrett estimat av den sanne effekten av log BNP på lykke.

Man kan videre inkludere region- og tidsfaste effekter i modellen, som i modell (2). Ved å justere for regionale faste effekter, fjerner vi all variasjon i lykkenivåer som kan forklares av region-spesifikke, tid-uavhengige faktorer. Dette er all konstant påvirkning fra for eksempel geografiske eller historiske forhold som er unike for hver region og som kan påvirke lykkenivået. Videre, ved å kontrollere for tidstrender, tar vi hensyn til generelle endringer i lykkenivåer over tid, uavhengig av region. Dette kan inkludere globale økonomiske trender eller teknologiske fremskritt som påvirker lykkenivåer over hele verden. Hensikten med disse faste effektene er å kontrollere for uobserverte variabler som kan variere mellom regioner og over tid, men som er konstante innen en region over tid (regionfaste effekter) eller som påvirker alle regioner likt over tid (tidsfaste effekter). Tilsvarende resonnement gjelder for siste modell (3), men hvor det er land-faste effekter istedenfor regioner. Når vi bruker en modell som inkluderer tid- og land-faste effekter så betyr det at vi bruker variasjon fra sammenhengen mellom lykke og log BNP som kommer fra endringer innad i et land over tid, som ikke er forklart av globale trender eller de andre forklaringsvariablene, til å estimere regresjonskoffisienten.

Likevel vil tidsvarierende forskjeller i land gjenstå. Hvis det fins variabler som vi ikke har kontrollert for som korrelerer både med utfallsvariabelen vår lifeladder og en forklaringsvariabel, for eksempel logGDP, vil vi ikke kunne tolke koffisienten til logGDP som den kausale effekten på lykke. Vi kan se for oss at det er en uobservert variabel som korrelerer med både BNP og lykke som varierer over tid i en region, som dermed vil være utelatt fra regresjonen. Et eksempel på en slik variabel kan være velferdspolitikk der økt BNP henger sammen med økte velferdsgoder som foreldrepermisjon, som igjen øker lykkenivåene i et land. Dette vil lede til forventningsskjeve estimater av effekten av log BNP på lykke. Dette vil dermed være et brudd på antagelsen om eksogenitet, der eksogenitet er definert som korrelasjon mellom forklaringsvariabelen og feilleddet.

Kommentar: Det er viktig at vi ber R tolke year som en kategorisk variabel ved hjelp av as.factor() kommandoen i regresjonen. Hvis ikke vil den tolkes som en nummerisk variabel og ikke inkludere dummy-variabler per år.

b) Diagnostiser en av regresjonsmodellene fra oppgave a).

Har her valgt å diagnostisere modell (2) med region- og tidsfaste effekter, men man kan så klart ha valgt å diagnostisere en annen modell.

Residualplott: Er det hetroskedastisitet?

residualer <- data.frame(residualer = fe_model_region$residuals,
                         predikert = fitted(fe_model_region))

ggplot(residualer, aes(x = predikert, y = residualer)) +
  geom_point()  +
  geom_hline(yintercept = 0, linetype = "dashed")+
  xlab("Predikert verdi")+
  ylab("Observert verdi")+
  ggtitle("Residualplott")+
  theme_classic()

Her ser vi at det ikke ser ut til å være heteroskedastisitet i residualene til verdier av predikert lifeladder.

Hva så med normalantagelsen om residualene?

fitted_normal <-
  tibble(x = seq(min(fe_model_region$residuals), max(fe_model_region$residuals), length.out = 200)) %>%
  mutate(y = dnorm(x, mean =0, sd = sd(fe_model_region$residuals)))

ggplot(residualer, aes(x = residualer, y = after_stat(density))) +
  geom_histogram()+
  geom_line(aes(x = x, y = y), data = fitted_normal, linetype = "dashed")+
  xlab("Residualer")+
  ylab("")+
  ggtitle("Histogram for observerte residualer")+
  theme_classic()

ggplot(residualer, aes(sample = residualer)) +
  stat_qq() +
  stat_qq_line() +
  xlab("Teoretiske kvantiler")+
  ylab("Observerte kvantiler")+
  ggtitle("QQ-plot")+
  theme_classic()

Fra histogrammet fremstår residualene normalfordelte, men med flere residualer sentrert nær null enn forventet. I QQ-plottet vil observasjonene ligge på den rette linjen hvis de er normalfordelte. Det at vi kan se en S-form på de plottede residualene betyr at residualene har tyngre haler enn det vi forventer dersom de er normalfordelte. Pga dette så vi skal være litt forsiktige med å stole på dekningsgraden til prediksjons- og konfidensintervaller.

Alt i alt fremstår det som regresjonsmodellen, modell (2), med tid- og regionsfaste effekter passer dataene våre relativt greit, og er en god modell for datasettet.

Oppgave 4: Last inn nytt datasett med en ny variabel og gjør flere analyser

Vi skal nå se på datasettet whr_ny. Dette datasettet inkluderer de balanserte WHR dataene fra whr_balansert og en ny variabel corruption. Denne variabelen er samlet inn av et FN organ (UNODC) og angir offisiell statistikk på antall korrupsjonssaker per 100,000 innbygger. Det vil si at vi nå har to mål på korrupsjon; oppfattet korrupsjon målt via spørreundersøkelser gjennomført av Gallup (percepcorruption) og offisiell statistikk rapportert til FN om antall korrupsjonssaker per 100,000 innbygger (corruption).

a) Hvilke sammenhenger finner vi mellom de to målene på korrupsjon? Illustrer med figur eller tabell. Hvilke inntrykk får dere av de to målene? Hva overrasker dere, og hva kan det bety?

ggplot(data = whr_ny, aes(x = percepcorruption, y = corruption)) +
  geom_point() +
  geom_smooth(method="lm", formula= y~x) +
  xlab("Perception of corruption") +
  ylab("Reported corruption cases per 100,000 inhabitants") +
  theme_classic()
Warning: Removed 669 rows containing non-finite values (`stat_smooth()`).
Warning: Removed 669 rows containing missing values (`geom_point()`).

whr_ny <- whr_ny %>%
  mutate(highlight_country = case_when(
    country %in% c("Norway", "Sweden", "Denmark") ~ country,
    TRUE ~ "Other"
  ))


ggplot(data = whr_ny, aes(x = percepcorruption, y = corruption, color = highlight_country)) +
  geom_point() +
  scale_color_manual(values = c("Sweden" = "yellow", "Denmark" = "red", "Other" = "gray")) +
  labs(color = "Country")+
  xlab("Perception of corruption")+
  ylab("Reported corruption cases per 100,000 inhabitants")+
  theme_classic()
Warning: Removed 669 rows containing missing values (`geom_point()`).

look <- whr_ny %>% 
  select(country, corruption, percepcorruption)

Hvis vi plotter de to korrupsjonsmålene ved hjelp av et scatterplot ser vi en overraskende negativ sammenheng mellom oppfattet korrupsjon og antall korrupsjonssaker per 100,000 innbygger.

Dette indikerer at offesiell statistikk ikke samsvarer med oppfattet korrupsjonsnivå. Når man bruker data bør man alltid tenke igjennom hvordan innsamlingmetoden vil kunne påvirke dataene. Noe som har gitt opphav til det presise statistikk begrepet “shit in, shit out”.

I dette tilfellet vil den litt rare sammenhengen kunne forklares med at noen land har interesse av å ikke rapportere sannferdig. Vi kan tenke oss at noen autoritære regimer er mindre villig til å rapportere inn korrupsjonssaker fordi det får dem til å fremstå dårlig, men det kan også tenkes at hva loven definerer som korrupsjon ikke er lik mellom landene. Derfor vil det i dette tilfellet trolig være bedre å bruke surveydata på korrupsjon enn offesiell statistikk som mål på faktisk korrupsjon.

Dette inntrykket blir forsterket når vi ser at det er overraskende nok Sverige og Danmark som har mange korrupsjonssaker, indikert i gult og rødt, to land som vi ofte ikke tenker på som veldig korrupte.

Refleksjon rundt offentlig statistikk på korrupsjon og lands styresett verdsettes. I dette tilfellet er nok surveydata mer troverdig enn offentlig statistikk, spesielt for regimer som ikke vil sannferdig rapportere slike hendelser.

b) Analyser sammenhengen mellom lykke (som avhengig variabel) og målene på korrupsjon, samt andre forklaringsvariabler ved bruk av regresjonsmodeller både med og uten faste effekter. Tolk resultatet.

# Gjør en regresjonsanalyse. Hva ser vi?
model1 <- lm(lifeladder ~  corruption + logGDP + socialsupport + freedom  + lifeexpectancy + generosity, data = whr_ny)
model2 <- lm(lifeladder ~  percepcorruption  + logGDP + socialsupport + freedom  + lifeexpectancy + generosity, data = whr_ny)
model3 <- lm(lifeladder ~  percepcorruption + corruption + logGDP + socialsupport + freedom  + lifeexpectancy + generosity, data = whr_ny)


model4 <- lm(lifeladder ~ corruption + logGDP + socialsupport + freedom  + lifeexpectancy + generosity + country + as.factor(year),
                        data = whr_ny)

model5 <- lm(lifeladder ~ percepcorruption + logGDP + socialsupport + freedom  + lifeexpectancy + generosity + country + as.factor(year),
                        data = whr_ny)

model6 <- lm(lifeladder ~ percepcorruption + corruption + logGDP + socialsupport + freedom  + lifeexpectancy + generosity + country + as.factor(year),
                        data = whr_ny)

stargazer(model1, model2, model3, model4, model5, model6, type="text", omit = c("country", "region", "year"))

=========================================================================================================================================================================
                                                                                     Dependent variable:                                                                 
                    -----------------------------------------------------------------------------------------------------------------------------------------------------
                                                                                         lifeladder                                                                      
                              (1)                      (2)                      (3)                      (4)                      (5)                      (6)           
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------
corruption                  0.003***                                           0.0001                   -0.001                                            -0.001         
                            (0.001)                                           (0.001)                  (0.001)                                           (0.001)         
                                                                                                                                                                         
percepcorruption                                    -0.974***                -1.641***                                         -1.009***                  -0.247         
                                                     (0.125)                  (0.168)                                           (0.221)                  (0.315)         
                                                                                                                                                                         
logGDP                      0.503***                 0.380***                  0.164*                  2.898***                 1.748***                 2.905***        
                            (0.084)                  (0.046)                  (0.091)                  (0.469)                  (0.185)                  (0.498)         
                                                                                                                                                                         
socialsupport               2.229***                 2.546***                 2.626***                 2.381***                 1.998***                 2.282***        
                            (0.498)                  (0.234)                  (0.467)                  (0.470)                  (0.325)                  (0.495)         
                                                                                                                                                                         
freedom                     2.255***                 1.330***                 0.938***                 0.761***                  0.366*                  0.756**         
                            (0.262)                  (0.181)                  (0.270)                  (0.291)                  (0.218)                  (0.300)         
                                                                                                                                                                         
lifeexpectancy              0.026**                  0.034***                 0.051***                  -0.054                 -0.101***                  -0.056         
                            (0.012)                  (0.006)                  (0.012)                  (0.052)                  (0.023)                  (0.053)         
                                                                                                                                                                         
generosity                  0.618***                 0.432***                  0.138                   0.550**                  0.542***                  0.514*         
                            (0.188)                  (0.127)                  (0.171)                  (0.265)                  (0.184)                  (0.271)         
                                                                                                                                                                         
Constant                   -4.602***                -2.666***                  -0.906                 -21.073***                -5.275**                -20.706***       
                            (0.492)                  (0.323)                  (0.571)                  (4.564)                  (2.193)                  (4.751)         
                                                                                                                                                                         
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Observations                  298                      797                      289                      298                      797                      289           
R2                           0.750                    0.778                    0.811                    0.960                    0.924                    0.959          
Adjusted R2                  0.745                    0.776                    0.806                    0.951                    0.914                    0.949          
Residual Std. Error     0.434 (df = 291)         0.495 (df = 790)         0.375 (df = 281)         0.191 (df = 238)         0.308 (df = 701)         0.193 (df = 229)    
F Statistic         145.683*** (df = 6; 291) 460.636*** (df = 6; 790) 172.081*** (df = 7; 281) 97.996*** (df = 59; 238) 89.648*** (df = 95; 701) 90.994*** (df = 59; 229)
=========================================================================================================================================================================
Note:                                                                                                                                         *p<0.1; **p<0.05; ***p<0.01

I modell (1) og modell (2) kjører vi en enkel (naiv?) regresjonsmodell hvor vi finner at antall korrupsjonssaker er positivt korrelert med lykke (modell (1)), mens oppfattet korrupsjon har en negativ assosiasjon med lykke (modell (2)). Begge statistisk signifikant forskjellig fra 0 på 1% nivå. Om vi hadde avsluttet analysen her, ville vi endt på en litt rar konklusjon. Hvis vi så kjører begge i samme modell (modell (3)), ser vi at det kun er oppfattet korrupsjon som gjenstår som statistisk signifikant, men antall korrupsjonssaker har fortsatt et positivt fortegn.

Dette inntrykket endrer seg noe når vi inkluderer land- og tidsfaste effekter for de tilsvarende modellene modell (4), modell (5) og modell (6). Her ser vi nemlig, til tross for manglende statistisk signifikans, at antall korrupsjonssaker får negativt fortegn.

Vi ser også at antall observasjoner i regresjonene varierer og går betydelig ned når vi inkluderer variabelen på antall korrupsjonsaker (pga. NAs). Derfor bør vi være forsiktige med å sammenligne på tvers av modeller fordi datagrunnlaget varierer.

Hvis vi sammenligner modell (1) og modell (4) så illustrerer det et viktig poeng om utelatt variabel bias, der det i modell (1) fremstår som en positiv statistisk signifikant sammenheng mellom korrupsjon og lykke om vi ikke tar hensyn til tidseffekter og, kanskje spesielt viktig, land-spesifikke effekter, som vi ser i modell (4).

Hvis det er slik at antall korrupsjonssaker rapportert av offisielle myndigheter egentlig er et slags mål på autoritære regimer (hvor lave verdier indikerer mer autoritært), vil en modell som holder land-spesifikke effekter fast faktisk kontrollere for det (gitt at det ikke endrer seg over tidsperioden vi ser på). Om man så har kontrollert for alt som er spesifikt for det landet, sitter man igjen med en koffisient på antall korrupsjonssaker som faktisk måler effekten på lykken, gitt at det ikke er andre utelatte variabler. Det at vi nå får et negativt fortegn indikerer at når et land, målt mot seg selv, har flere korrupsjonssaker så er det assisoiert med lavere lykke. Det er riktignok ikke statistisk signifikant, noe som kan være pga. det lave antallet observasjoner, høy variasjon, eller at det faktisk ikke er så sterk sammenheng mellom lykke og korrupsjonssaker.

Fra modell (6) ser vi at begge koffisientene til korrupsjonsmål er negative slik som vi ville antatt om sammenhengen mellom korrupsjon og lykke. Riktignok forsvinner signifikansen til begge korrupsjonsvariablene. Dette kan være pga. for stor variasjon i det begrensede datasettet hvor begge variablene eksisterer.