3.6 Relevante R-kommandoer

Under følger en liste over hvilke oppgaver du skal klare i R fra denne modulen. Vår policy er at R-kommandoene under er tilstrekkelige for å løse oppgavene i datalabber og eksamen 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

Di 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.

\(t\)-tester

Du må kunne gjøre de ulike \(t\)-testene ved hjelp av t.test()-funksjonen. vi bruker eksempeldatasettet testdata.xls som illustrasjon, som inneholder de to kolonnene X1 og X2. I bunn og grunn handler det om å hente ut de relevante kolonnene i datasettet som vektorer, og så å sette argumentene i t.test()-funksjonen riktig.

# Tosidig, ett-utvalgs t-test for om forventningen til X1 er lik 100:
t.test(testdata$X1, mu = 100)

# Samme testen, men ensidig med alternativhypotese om at mu > 100
t.test(testdata$X1, mu = 100, alternative = "greater")

# For å gjøre ensidig test den andre veien bytter du ut "greater" med "less".

# Vanlig to-utvalgs t-test for om forventningen til X1 og X2 er like:
t.test(testdata$X1, testdata$X2)

# Vanlig to-utvalgs t-test der vi antar at variansen i de to populasjonene er like:
t.test(testdata$X1, testdata$X2, var.equal = TRUE)

# Parret to-utvalgs t-test, som bare gir mening dersom de to vektorene med
# observasjoner er like lange, og dersom observasjonene er matchet opp slik at
# første element i den første vektoren skal matches mot første element i den
# andre vektoren etc.:
t.test(testdata$X1, testdata$X2, paired = TRUE)

Vi kan bytte ut alternativhypotesen til ensidige tester over alt ved å sette alternative = til ønskelig verdi.

Ensidige to-utvalgs tester: Hvilken vei tester vi?

Når du gjør en ensidig to-utvalgs \(t\)-test er det viktig å forstå hva alternative = "greater" og alternative = "less" faktisk betyr. Funksjonen t.test() ser alltid på differansen mellom første og andre argument. Altså:

  • t.test(x, y, alternative = "greater") tester \(H_1: \mu_x - \mu_y > 0\), dvs. at gjennomsnittet i x er høyere enn i y.
  • t.test(x, y, alternative = "less") tester \(H_1: \mu_x - \mu_y < 0\), dvs. at gjennomsnittet i x er lavere enn i y.

Det betyr at du kan snu på resultatet ved å bytte rekkefølgen på argumentene. Følgende to tester er ekvivalente:

t.test(testdata$X1, testdata$X2, alternative = "greater")
t.test(testdata$X2, testdata$X1, alternative = "less")

Tips: Sjekk alltid R-utskriften! Den skriver eksplisitt ut hva alternativhypotesen er, for eksempel: "alternative hypothesis: true difference in means is greater than 0". Da ser du at testen samsvarer med det du ønsker å teste.

Bruke formel-syntaks med ~

I mange praktiske situasjoner har vi et datasett der en kategorisk variabel angir gruppemedlemskap. Da kan vi bruke formelsyntaksen variabel ~ gruppe i stedet for å filtrere ut to separate vektorer. La oss si at vi har et datasett df med en numerisk variabel verdi og en kategorisk variabel gruppe med nivåene "A" og "B":

t.test(verdi ~ gruppe, data = df)

Viktig: Når du bruker ~-syntaksen, bestemmes rekkefølgen av nivårekkefølgen til den kategoriske variabelen (alfabetisk som standard). Det betyr at R ser på differansen mellom det første nivået og det andre nivået. Dersom nivåene er "No" og "Yes", vil R teste differansen No - Yes, fordi "No" kommer før "Yes" alfabetisk. Utskriften fra t.test() vil alltid vise dette eksplisitt, for eksempel: "alternative hypothesis: true difference in means between group No and group Yes is less than 0". For ensidige tester må du derfor tilpasse alternative-argumentet til denne rekkefølgen: Dersom du ønsker å teste om gruppe "Yes" har høyere gjennomsnitt enn gruppe "No", tester du om differansen No - Yes er negativ, altså bruker du alternative = "less". Tilsvarende, hvis du tror at "No" har høyest gjennomsnitt, bruker du alternative = "greater".

Ulike utvalgsstørrelser

Videre kan vi tenke oss tilfeller der observasjonsvektorene ikke er like lange, og for eksempel kommer i en excel-fil slik som denne der den ene kolonnen inneholder flere tall enn den andre. Det gjør ingenting. Vi kan lese inn datasettet på vanlig måte, og legger merke til at de “tomme” plassene i tabellen blir fylt med NA. Vi kan kjøre \(t\)-testene på vanlig måte.

testdata2 <- readxl::read_excel("ulik-lengde.xlsx")
testdata2
t.test(testdata2$X1, testdata2$X2)
## # A tibble: 15 × 2
##       X1    X2
##    <dbl> <dbl>
##  1     1     5
##  2     2     6
##  3     3     7
##  4     4     8
##  5     5     9
##  6     6    10
##  7     7    11
##  8     8    12
##  9     9    13
## 10    10    14
## 11    11    NA
## 12    12    NA
## 13    13    NA
## 14    14    NA
## 15    15    NA
## 
##  Welch Two Sample t-test
## 
## data:  testdata2$X1 and testdata2$X2
## t = -1, df = 22.975, p-value = 0.3277
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4.603173  1.603173
## sample estimates:
## mean of x mean of y 
##       8.0       9.5

Varianstester

Varianstest for ett uvalg er ikke bygget inn i R som en egen funksjon. Da må vi gjøre den “semi-automatisk” og regne ut testobservatoren ved å bruke formelen direkte. Her tester vi om variansen til X2-kolonnen i testdatasettet er lik 10:

# Varianstest
s <- sd(testdata$X2)
sigma_0 <- 10
n <- length(testdata$X2)

# Testobservatoren
T <- (n-1)*s^2/sigma_0^2

# Kritiske verdier
L <- qchisq(0.025, df = n - 1)
U <- qchisq(0.975, df = n - 1)

Vi må så sammenligne verdien av testobservatoren T med de kritiske verdiene L og U. Legg merke til at vi har brukt en funksjon i R (qchisq()) for å finne de aktuelle kvantilene i kjikvadratfordelingen.

For to-utvalgs varianstest kan vi bruke funksjonen var.test():

var.test(testdata$X1, testdata$X2)

Vi kan også her gjøre ensidige tester ved å sette alternative =-argumentet til enten "less" eller "greater".

Tester for andeler

Både for ett-utvalgs og to-utvalgs test for andeler kan vi fylle inn tall direkte i testobservatoren på denne måten (sett inn for de nødvendige verdiene):

# Ett-utvalgs test for andeler
p_hatt <- mean(testdata$A1)
p_null <- 0.5
n      <- length(testdata$A1) 

Z <- (p_hatt - p_null)/sqrt(p_null*(1-p_null)/n)   

# To-utvalgs test for andeler
p1_hatt <- mean(testdata$A1)
p2_hatt <- mean(testdata$A2)

n1 <- length(testdata$X1)
n2 <- length(testdata$X2)

p_hatt <- mean(c(testdata$A1, testdata$A2))        # <- Regner ut felles p_hatt

Z <- (p1_hatt - p2_hatt)/sqrt((1/n1 + 1/n2)*p_hatt*(1-p_hatt))

Testobservatorene må da sammenlingnes med relevant kvantil i normalfordelingen, som du enten finner i tabell i læreboken eller ved å bruke funksjonen qnorm() med ønsket kvantil som argument.

Kjikvadrattester for goodness-of-fit

Hvis vi skal teste for goodness-of-fit for en fordeling kan vi raskt regne ut testobservatoren direkte. Vi bruker eksempelet fra forelesningsvideoen:

p0 <- c(0.45, 0.40, 0.15)  # Fordeling under H0
f  <- c(102, 82, 16)       # Observerte frekvenser
e  <- p0*sum(f)            # Forv. frekv. under H0

T  <- sum((f - e)^2/e)     # Testobservator

c  <- qchisq(.95, df = 2)  # Kritisk verdi

# Forkast H0?
T > c

Alternativt så kan vi bruke chisq.test()-funksjonen direkte. Pass på å bruke riktige argumentnavn (x = og p =) når du skal gjøre goodness-of-fit:

# Eventuelt direkte
chisq.test(x = f, p = p0)

Dersom du skal teste for uavhengighet så kan du anta at datasettet kommer som en ferdig kontingenstabell, med den ene gruppeinndelingen langs kolonnene, og den andre gruppeinndelingen langs radene. En slik tabell kan du sende rett inn i chisq.test()-funksjonen for å gjøre en test for uavhengighet:

# Test for uavhengighet
immigration <-
  read_xls("immigration-wide.xls",
           range = "B2:E12",
           col_names = c("many", "some", "few", "none"))
chisq.test(immigration)

Kjikvadrattester for uavhengighet

Som eksempel kan vi ta oppgaven fra Dataøving 2. Vi starter med å lese inn dataene:

# les inn data
library(readxl)
violence <- read_excel("violence.xlsx")

Vi bruker select() til å velge de to variablene vi er interessert i og funksjonen table() til å lage en krysstabell:

violence_redusert <-
  violence %>%
  select(violent_treatment, experienced_violence)
krysstabell <- table(violence_redusert)
krysstabell
##                  experienced_violence
## violent_treatment Less Violent Violent
##      Less Violent          114       9
##      Violent                33      93

Det ser ut til at det er en klar sammenheng mellom faktisk og opplevd voldelighet ved at de fleste forsøkspersonene havner på diagonalen i krysstabellen over.

Nå skal vi utføre en test av nullhypotesen \(H_0\): faktisk og opplevd voldelighet er uavhengige kjennetegn, mot den alternative hypotesen \(H_1:\) at de er avhengige kjennetegn. Dette er en chikvadrattest som kan utføres i R på følgende måte:

chisq.test(krysstabell)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  krysstabell
## X-squared = 111.06, df = 1, p-value < 2.2e-16

og vi ser at vi får en soleklar forkastning. Dette er også logisk gitt hva vi faktisk tester her.

Kritiske verdier i R

Som vi har sett over baserer tester i R seg stort sett på en eller annen datainput. Som output får du både testobservator og p-verdi. Hvis du allikevel av en eller annen grunn (f.eks. på eksamen) trenger å finne en kritisk verdi i en fordeling, så kan det være nyttig å vite hvordan man gjør dette i R.

De kritiske verdiene finner vi ved å bruke kvantilfunksjonene til den respektive fordelingen. Disse funksjonene har alltid et argument p og har som output det tallet som har sannsynlighetsmasse p til venstre for seg. Du må derfor velge p ut fra signifikansnivå og om det er en ensidig/tosidig test.

Kritiske verdier i T-fordelingen

Dersom du f.eks skal utføre en t-test for \(H_0: \mu = \mu_0\) mot \(H_1: \mu < \mu_0\), med signifikansnivå \(\alpha = 0.05\) og har \(40\) observasjoner er kritisk verdi:

qt(0.05, df = 40 - 1)
## [1] -1.684875

og du forkaster da dersom din testobservator \(T\) er mindre enn dette tallet.

Hadde alternativ hypotesen vært \(H_1: \mu > \mu_0\) ville kritisk verdi vært

qt(1 - 0.05, df = 40 - 1)
## [1] 1.684875

og du forkaster dersom din testobservator \(T\) er større enn dette tallet.

For en to-sidig test (\(H_1: \mu\neq \mu_0\)) forkaster du dersom \(T\) faller uten for intervallet:

L <- qt(0.05/2, df = 40 - 1)
U <- qt(1 - 0.05/2, df = 40 - 1)
c(U,L)
## [1]  2.022691 -2.022691

Kritiske verdier i standardnormalfordelingen

Helt likt som for t-fordelingen bare at det ikke er noen frihetsgrader involvert:

# Ensidig mu < mu0, alpha = 0.05
qnorm(0.05)
## [1] -1.644854
# Ensidig mu > mu0, alpha = 0.05
qnorm(1 - 0.05)
## [1] 1.644854
# Tosidig
L <- qnorm(0.05/2)
U <- qnorm(1 - 0.05/2)
c(U,L)
## [1]  1.959964 -1.959964

Alternativt er det bare å pugge tallene \(1.64\) og \(1.96\).

Kritiske verdier i \(\chi^2\)-fordelingen

\(\chi^2\)-fordelingen brukes i to ulike sammenhenger i MET4: varianstester og kjikvadrattester (goodness-of-fit og uavhengighet). For kjikvadrattestene er det kun den øvre halen som er relevant, fordi en stor testobservator betyr at de observerte frekvensene er langt fra de forventede. For varianstester kan vi derimot trenge nedre hale, øvre hale eller begge, avhengig av alternativhypotesen.

Varianstester:

Dersom du tester \(H_0: \sigma^2 = \sigma_0^2\) mot \(H_1: \sigma^2 < \sigma_0^2\) (variansen er mindre enn påstått), forkaster du i nedre hale:

# Kritisk verdi, 5% nivå, ensidig test, nedre hale, 30 frihetsgrader
qchisq(0.05, df = 30)
## [1] 18.49266

og du forkaster dersom testobservatoren er mindre enn dette tallet.

Dersom du tester mot \(H_1: \sigma^2 > \sigma_0^2\) (variansen er større enn påstått), forkaster du i øvre hale:

# Kritisk verdi, 5% nivå, ensidig test, øvre hale, 30 frihetsgrader
qchisq(1 - 0.05, df = 30)
## [1] 43.77297

og du forkaster dersom testobservatoren er større enn dette tallet.

For en tosidig test (\(H_1: \sigma^2 \neq \sigma_0^2\)) trenger du begge grensene:

# Kritisk verdi, 5% nivå, tosidig test, 30 frihetsgrader
L <- qchisq(0.05/2, df = 30)
U <- qchisq(1 - 0.05/2, df = 30)
c(L, U)
## [1] 16.79077 46.97924

og du forkaster nullhypotesen dersom testobservatoren havner utenfor dette intervallet.

Kjikvadrattester (goodness-of-fit og uavhengighet):

Her forkaster vi alltid i øvre hale, og trenger kun:

# Kritisk verdi, 5% nivå, k = 3 kategorier (k - 1 frihetsgrader)
qchisq(1 - 0.05, df = 2)
## [1] 5.991465

og du forkaster dersom testobservatoren er større enn dette tallet.

Kritiske verdier i F-fordelingen

F-fordelingen trenger vi når vi skal teste om to varianser er like. Skal du gjøre denne testen manuelt er det lurt å sette den største variansen i telleren av testobservatoren

\[F = \frac{S_{1}^2}{S_{2}^2}.\]

Uavhengig om det er en ensidig eller tosidig test, trenger vi da bare å sjekke om testobservatoren overstiger den kritiske verdien i den øvre halen i F-fordelingen.

Har vi brukt \(30\) observasjoner til å regne ut \(S_{1}^2\) og \(25\) observasjoner til å regne ut \(S_{1}^2\) blir kritisk verdi i en to-sidig test der \(H_1: \sigma_{1}^2 \neq \sigma_{2}^2\) :

# Kritisk verdi, 5% nivå, tosidig test

qf(1 - 0.05/2, df1 = 30 - 1, df2 = 25 - 1)
## [1] 2.217443

Er det en ensidig test hvor \(H_1: \sigma_{1}^2 > \sigma_{2}^2\) blir kritisk verdi:

# Kritisk verdi, 5% nivå, ensidig test

qf(1 - 0.05, df1 = 30 - 1, df2 = 25 - 1)
## [1] 1.945259

Husk: Du kan definere hva som er \(S_{1}^2\) og \(S_{2}^2\), og med strategien over definerer du alltid \(S_{1}^2\) til være den av variansene som er størst.