4.5 Relevante R-kommandoer

Under følger en liste over hvilke oppgaver du skal klare i R fra denne modulen. Vår policy fra og med vårsemesteret 2022 er at R-kommandoene under er tilstrekkelige for å løse oppgavene i datalabber og hjemmeeksamen i MET4. Det er med andre ord ikke nødvendig å lære seg teknikker utover det som er listet opp eksplisitt i listen under. Eventuelle nye teknikker som trengs for å løse en bestemt oppgave vil bli oppgitt og forklart dersom det er nødvendig. Det antas i tillegg at du kan den grunnleggende R-syntaksen som er dekket under Introduksjon til R.

Antakelser om datasett

Du kan gjøre de samme antakelsene om datasettet som i den tilsvarende oversikten i forrige modul: Datasettene er inneholdt i Excel-filer (.xls eller xslx) eller .csv-filer, som kolonner av variabler, med variabelnavn i første rad.

Tilpasse en lineær regresjonsmodell

Her er et lekedatasett med tre forklaringsvariabler x1, x2 og x3, samt en responsvariabel y: data-reg.xsls.

# Last inn data
library(readxl)
df <- read_excel("data-reg.xlsx")

# Syntaksen for å tilpasse en lineær regresjonsmodell med responsvariabel y og
# forklaringsvariabler x1, x2, x3 i datasettet df er som følger:
reg1 <- lm(y ~ x1 + x2 + x3, data = df)

# Dersom vi vil ha med interaksjonsledded mellom x1 og x3 kan vi skrive
# kommandoen under. Da vil også x1 og x3 automatisk bli med som variabler i
# tillegg til interaksjonsledde x1*x3.
reg2 <- lm(y ~ x1 + x2*x3, data = df)

# Log transformasjoner kan vi skrive rett inn i formelen
reg3 <- lm(log(y) ~ x1 + log(x2) + x3, data = df)

# Andre transformasjoner, som for eksempel dersom vi ønsker å ta med
# andregradsleddet x1^2, er greiest å få til dersom vi først lager en ny kolonne
# med kvadratene til x1, og så lager til regresjonen på vanlig måte:
df_ny <- df %>% mutate(x1_kvadrat = x1^2)
reg4 <- lm(y ~ x1 + x1_kvadrat+ x2 + x3, data = df_ny)

# Vi kan skrive ut et sammendrag for regresjonsmodellen ved hjelp av summary():
summary(reg1)

# Pakken stargazer inneholder en funksjon for å lage pene regresjonstabeller:
library(stargazer)
stargazer(reg1, reg2, type = "text")

# Man kan endre type-argumentet til "html", og sette out-argumentet til et
# filnavn hvis du vil skrive ut tabellen til en fil.

Prediksjoner, med konfidens- og prediksjonsintervall

La oss si at vi ønsker å predikere verdien av responsvariabelen for følgende verdier av forklaringsvariablene:

x1	x2	x3
10	95	1
20	96	1
30	100	0

Da må vi først samle verdiene av forklaringsvariablene i en data frame; det kan vi enten gjøre direkte i R:

ny_data <- data.frame(x1 = c(10, 20, 30), x2 = c(95, 96, 100), x3 = c(1,1,0))

Et alternativ kan være å lage til et regneark i Excel, og så lese det inn ved hjelp av read_excel(). Det er viktig at variabelnavnene i det nye datasettet er eksakt de samme som forklaringsvariablene i datasettet du brukte til å estimere regresjonsmodellen! Vi kan lage prediksjon med 95% konfidensintervall ved hjelp av predict():

prediksjoner <- predict(reg1, 
                        newdata = ny_data, 
                        interval = "confidence", 
                        level = 0.95)

Resultatet blir en ny data frame med kolonnene fit (prediksjoner), samt lwr og upr, som inneholder nedre og øvre grense i konfidensintervallet henholdsvis. Du kan endre til prediksjonsintervaller ved å bytte ut "confidence" med "prediction" (eventuelt "none" hvis du ikke ønsker noe intervall).

Lage residualplott

Vi må kunne lage følgende residualplott i R, her demonstrert med regresjonsmodellen reg1 som vi laget over:

# Henter ut residualer og predikerte verdier fra modellen, og setter de inn som 
# kolonner i en dataframe:
residualer <- data.frame(residualer = reg1$residuals,
                         predikert = reg1$fitted.values)

# Spredningsplott, residualer mot predikerte verdier
ggplot(residualer, aes(x = predikert, y = residualer)) +
  geom_point()

# Autokorrelasjonsplott
acf(reg1$residuals)

# Histogram for å sjekke normalantakelsen
ggplot(residualer, aes(x = residualer)) +
  geom_histogram()

# QQ-plott, observasjonene ligger langs en rett linje hvis de er normalfordelte
ggplot(residualer, aes(sample = residualer)) +
  stat_qq() +
  stat_qq_line()

# Regner ut og plotter Cooks avstand:
infl <- influence.measures(reg1)

# Lager enkelt plott av Cooks avstand
cooks_avstand <- data.frame(obs = 1:nrow(df),
                            cook = infl$infmat[, "cook.d"],
                            cook_infl = infl$is.inf[, "cook.d"])

ggplot(cooks_avstand, aes(x = obs, y = cook)) +
  geom_col()

Alle plottene over kan justeres og pyntes, for eksempel ved å gjøre endringer som i R-videoen om plotting. Se også forelesningsscriptet for flere tips om presentasjon av residualplott.