library(dplyr)
library(ggplot2)
load("met4_h22.Rdata")

Denne hjemmeksamen har ikke et ‘korrekt’ løsningforslag ettersom flere problemstillinger kan formuleres og undersøkes, f.eks:

Nøkkelen til en god besvarelse er at det skrives en konsis analyse med en rød tråd og med bruk av et begrenset antall metoder. Å begrense oppgaven til delmenger av dataene kan være fornuftig, men en får også færre data å basere seg på så her er det en fin balansegang.

For å gi input til flest mulig vil vi under gjøre noen generelle betraktninger om datasettet som kan være relevante for flere av problemstillingene over.

Et viktig poeng omtrent uansett hvilken problemstilling en går for er at en må ta inn over seg at variabelen loss er en andel. Den kan rett og slett ikke være normalfordelt siden den er begrenset til å ligge mellom 0 og 1. Dette ser vi lett dersom vi lager et histogram over variabelen:

wastedata %>% 
  select(loss) %>% 
  ggplot(aes(x = loss)) + geom_histogram()

eller ser på en kvantitativ oppsummering av denne variabelen:

wastedata %>% 
  select(loss) %>% summary()
##       loss          
##  Min.   :0.0000014  
##  1st Qu.:0.0146000  
##  Median :0.0301000  
##  Mean   :0.0567311  
##  3rd Qu.:0.0607500  
##  Max.   :0.9862000

Siden minste observerte verdi er større enn 0 kan vi logtransformere denne variabelen og se om vi får noe som ser mer normalfordelt ut:

wastedata %>% 
  select(loss) %>% 
  mutate(log_loss = log(loss)) %>% 
  ggplot(aes(x = log_loss)) + geom_histogram()

Det er med andre fornuftig å bruke log(loss) som responsvariabel når man utfører forskjellige regresjonsanalyser.

Det vil generelt være interessant å undersøke hvordan mattap varierer med år, land og korntype. En forenkling som er lur å gjøre innledningsvis er at en ser på en av disse variablene om gangen ved hjelp av group_by() og summarize() som forklart i forberedelsesinformasjonen. Det vi hovedsak tar med oss plottene under er at korntapet varierer i stor grad mellom land og korntype, men (desverre) ikke så mye i tid selv om vi ser en liten nedgang i slutten av perioden.

En figur med et stolpediagram av gjennomsnittlig tap per land kan man f.eks lage ved koden:

wastedata %>% 
  group_by(country) %>% 
  summarize(loss = mean(loss)) %>% 
  ggplot(aes(x = reorder(country,-loss), y = loss)) +  # 'reorder' for å få fra størst til minst 
  geom_bar(stat = "identity") + 
  theme(axis.text.x = element_text(angle=90, vjust=.5, hjust=1)) + # kan være lurt å rotere navnene på x-aksen
  xlab("land")

evt. et boksplott for å se litt på fordelingen i hvert land:

wastedata %>% 
  ggplot(aes(x = reorder(country,-loss), y = loss)) +  # 'reorder' for å få fra størst til minst 
  geom_boxplot() + 
  theme(axis.text.x = element_text(angle=90, vjust=.5, hjust=1)) + # kan være lurt å rotere navnene på x-aksen
  xlab("land")

Her ser vi at det stor variasjon mellom land. Tilsvarende for gjennomsnittlig tap per korntype:

wastedata %>% 
  group_by(crop) %>% 
  summarize(loss = mean(loss)) %>% 
  ggplot(aes(x = reorder(crop,-loss), y = loss)) +  
  geom_bar(stat = "identity") + 
  theme(axis.text.x = element_text(angle=90, vjust=.5, hjust=1)) +
  xlab("korntype")

Vi kan gjøre det samme for år men da er det naturlig å ikke ordne årstallene etter størrelse på tap, men la dem stå i kronologisk rekkefølge:

wastedata %>% 
  group_by(year) %>% 
  summarize(loss = mean(loss)) %>% 
  ggplot(aes(x = year, y = loss)) +  
  geom_bar(stat = "identity") + 
  theme(axis.text.x = element_text(angle=90, vjust=.5, hjust=1)) +
  xlab("årstall")

Før en går i gang med en modell er det alltid lurt å se litt på histogram, spredningsplott og korrelasjoner mellom variablene. Dette kan gjøres enkeltvis/parvis, men her tar vi en liten snarvei ved hjelp av ggpairs() funksjonen i pakken GGally (ikke pensum, men gjør denne rapporten mye kortere):

library(GGally)
wastedata %>% 
  filter(crop == "wheat") %>% 
  select(-c(country, year, crop)) %>%
ggpairs()

Her har vi hentet ut tapet av “wheat” og sett på korrelasjon, spredningsplott, og fordelingene per variabel. Vi har allerede diskutert at det kan være lurt å log-transformere “loss”, men vi ser ut fra diagonalplottene i figuren over at også at biofuel, gas, GDP og coal med fordel kan log-transformeres for å få bedre kontroll på uteliggere. Siden disse verdiene tar negative verdier plusser vi på en stor nok konstant (i dette tilfellet holder det med 1) før vi tar logaritmen. For grupper som har fokusert på visse forklaringsvariabler er det viktig å se på korrelasjonen mellom de inkluderte og utelatte variablene (“omitted variable bias”).

Transformasjoner kan vi gjøre direkte i modellbyggingen. Skal vi f.eks.tilpasse en paneldatamodell for log-tap av hvete kan vi gjøre følgende:

library(plm)
waste_panel <- wastedata %>% 
  filter(crop == "wheat") %>%
  pdata.frame(index = c("country", "year")) 


waste_model <- plm(log(loss) ~ year + temperature + rain + log(biofuel + 1) + log(gas + 1) + 
                     log(GDPcontribution + 1) + log(coal + 1) + lpi + input + investment + energy,
                   data = waste_panel,
                   model = "within")
summary(waste_model)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = log(loss) ~ year + temperature + rain + log(biofuel + 
##     1) + log(gas + 1) + log(GDPcontribution + 1) + log(coal + 
##     1) + lpi + input + investment + energy, data = waste_panel, 
##     model = "within")
## 
## Unbalanced Panel: n = 41, T = 8-24, N = 828
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -1.8014990 -0.1697518  0.0084393  0.1644840  2.3886066 
## 
## Coefficients:
##                            Estimate Std. Error t-value Pr(>|t|)  
## year1992                  0.1813687  0.1095047  1.6563  0.09808 .
## year1993                  0.1700970  0.1104484  1.5401  0.12397  
## year1994                  0.1842877  0.1101367  1.6733  0.09469 .
## year1995                  0.1130710  0.1098817  1.0290  0.30380  
## year1996                  0.1820450  0.1097655  1.6585  0.09763 .
## year1997                  0.1308721  0.1097766  1.1922  0.23357  
## year1998                  0.1122052  0.1100863  1.0192  0.30841  
## year1999                  0.1406854  0.1105727  1.2723  0.20365  
## year2000                  0.1319134  0.1125019  1.1725  0.24135  
## year2001                  0.1513479  0.1122800  1.3480  0.17808  
## year2002                  0.1971071  0.1120404  1.7593  0.07894 .
## year2003                  0.1949185  0.1134080  1.7187  0.08607 .
## year2004                  0.1632495  0.1184848  1.3778  0.16867  
## year2005                  0.1765454  0.1244003  1.4192  0.15626  
## year2006                  0.2720995  0.1502541  1.8109  0.07055 .
## year2007                  0.2511974  0.1850902  1.3572  0.17514  
## year2008                  0.2545277  0.2383510  1.0679  0.28592  
## year2009                  0.3309897  0.1848099  1.7910  0.07370 .
## year2010                  0.3828873  0.2079270  1.8415  0.06595 .
## year2011                  0.3244443  0.2341316  1.3857  0.16624  
## year2012                  0.3404628  0.2356715  1.4446  0.14897  
## year2013                  0.2573675  0.2180007  1.1806  0.23814  
## year2014                  0.1990566  0.2202830  0.9036  0.36647  
## temperature               0.0309635  0.0323005  0.9586  0.33806  
## rain                      0.0619862  0.0368262  1.6832  0.09275 .
## log(biofuel + 1)          0.0537030  0.1212290  0.4430  0.65790  
## log(gas + 1)             -0.1805417  0.0998424 -1.8083  0.07096 .
## log(GDPcontribution + 1) -0.0053252  0.0691549 -0.0770  0.93864  
## log(coal + 1)             0.0138924  0.0581781  0.2388  0.81133  
## lpi                      -0.0819442  0.0370520 -2.2116  0.02729 *
## input                    -0.0212148  0.0737164 -0.2878  0.77359  
## investment                0.0042336  0.0472191  0.0897  0.92858  
## energy                    0.0537595  0.0338385  1.5887  0.11255  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    160.15
## Residual Sum of Squares: 153.72
## R-Squared:      0.040145
## Adj. R-Squared: -0.052785
## F-statistic: 0.955627 on 33 and 754 DF, p-value: 0.54058

Her ser vi at for hvete så finner vi ingen signifikante forklaringsvariabler (Se på F-testen), etter vi har justert for faste land- og årseffekter. Da gir det heller ingen mening å begynne å tolke de forskjellige koeffisientene.

Resultatene vil selvsagt variere om en ser på spesifikke land, andre korntyper, eller andre tidsperioder og vi vil ikke undersøke alle varianter her.

Vi tar selvkritikk når det gjelder kursets dekning av diagnostisering av paneldatamodeller og vi har derfor vært ganske tilgivende på spørsmål 4 i oppgaven. Men fremgangsmåten er svært lik som ved vanlig regresjon. Vi kan hente ut residualene på samme måte:

residualer <- data.frame(residualer = waste_model$residuals)

og lage histogram og QQ-plott på vanlig måte:

ggplot(residualer, aes(x = residualer)) +
  geom_histogram()

ggplot(residualer, aes(sample = residualer)) +
  stat_qq() +
  stat_qq_line()

Vi ser at det er en del ekstreme observasjoner som gjør at normalantagelsen ikke holder. Husk: Diagnostiseringen trenger ikke være en “obduksjon” av modellen. Når en ser slike figurer er det naturlig å prøve finne og evt. utelukke de ekstreme observasjonene av loss for så å kjøre modellen på nytt. Tilsvarende figurer uten log-transformering ser enda verre ut og det samme gjelder da.

Det er helt greit at det ikke er funnet noen signifikante forklaringsvariabler eller at resultatet ikke gir en klar konklusjon på problemstillingen som er formulert, noe som var tilfelle for flere grupper på denne eksamen. Da kan det være lurt å diskutere hvorfor i stedet for å tvinge gjennom tolkninger som ikke støttes opp av analysen.