3.6 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

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 elemenet i den
# andre vektoret 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.

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. Alternativt kan vi slå opp i tabell.

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å en hjemmeeksamen) trenger å finne en kritisk verdi i en fordeling, så kan det være nyttig å vite hvordan man gjør dette i R i stedet for å slå opp i en eller annen tabell i en bok.

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

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

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

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

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.