Hvad er normalt? :: Draft

  Claus Thorn Ekstrøm    18. mar 2018    5 min

Troværdigheden af resultaterne fra en statistisk analyse står og falder med om modelantagelserne er opfyldte. Er de underliggende antagelser ikke opfyldte, kan vi ikke være sikre på estimater, konfidensintervaller og p værdier fra analysen. Et spørgsmål, der ofte bliver stillet i forbindelse med statistiske konsultationer er, om data er normalfordelte. Men hvad er det, der skal være normalfordelt og hvor kritisk er det egentlig, hvis data ikke opfylder normalfordelingsantagelserne?

De fleste har hørt, at hvis betingelserne omkring normalfordelingen ikke passer, så skal man gøre “noget andet”, hvilket oftest består i at lave en eller anden form for simpel ikke-parametrisk statistik.

Ulempen ved disse ikke-parametriske metoder er, at de primært beskæftiger sig med test af hypoteser, og at fokus derfor flyttes fra at diskutere og fortolke effektestimater til udelukkende at se på p værdier. Med andre ord bliver det kun muligt at vurdere, om et fund er statistisk signifikant og ikke direkte muligt at vurdere, om et eventuelt fund er videnskabligt relevant. Flere af de mest almindelige ikke-parametriske metoder har desuden den ulempe, at de ikke altid så let kan modellere effekter justeret for andre forklarende variable samtidig.1 Der er derfor gode grunde til at holde sig til de parametriske modeller, og normelfordelingsantagelserne er heldigvis så robuste, at man i mange tilfælde kan få fornuftige og tilfredsstillende resultater ud af lineære modeller også selvom normalfordelingsantagelserne ikke direkte er opfyldte.

Hvad er det, der skal være normalfordelt?

Lineære modeller (som omfatter specialtilfælde som t test, lineær regression, multipel lineær regression, ANOVA, ANCOVA) er en utrolig flektisbel klasse af modeller til statistisk analyse, og her illustreres betydningen af normalfordelingsantagelserne gennem en simpel lineær model med bare en enkelt forklarende variabel, \(x\). Den lineære model er givet ved formlen

\[y_i = \alpha + \beta\cdot x_i + \varepsilon_i,\]

hvor \(y_i\) og \(x_i\) er henholdsvis udfaldet og forklarende variabel for den \(i\)’te observation, \(\alpha\) er den gennemsnitlige værdi, når \(x=0\), \(\beta\) er den forventede ændring i gennemsnittet, når \(x\) stiger med 1 og \(\varepsilon\) er den “støj”, der er rundt om regressionslinjen. For lineære modeller antages det, at støjen er normalfordelt (med middelværdi 0 og konstant varians), \(\varepsilon \sim N(0, \sigma^2)\).

Lad os starte med at slå modellens antagelser fast: det er ikke data i sig selv, der skal være normalfordelte, men støjen. Det er \(\varepsilon\), der skal være normalfordelt, og det giver derfor ikke mening at lave histogrammer, qq-plots, eller normalfordelingstests for hverken \(x\)’erne i datasættet eller \(y\)’erne i datasættet (se figur ??).2 Modelkontrollen skal baseres på residualerne, \(\varepsilon\).

Figur 1: Tre eksempler på sammenhænge, der alle opfylder normalfordelingsantagelserne. Over hvert plot er vist et histogram for \(x\)’erne og ved siden er et histogram over \(y\)’erne. Normalfordelingsantagelserne er opfyldte for alle tre plots, men alligevel gælder der, at ingen af \(x\)’erne eller \(y\)’erne er normalfordelte.

Check af normalfordelingsantagelser

Der er forskellige … til at checke

\[\hat\varepsilon_i = y_i - (\hat\alpha + \hat\beta \cdot x_i)\]

qq plot

residualplot

Varianshomogenitet

Er robust

CLT - store stikprøver

library("nortest")

N <- c(20, 50, 100, 200)

runsim <- function(N) {

B <- 1000

res <- sapply(1:B, function(i) {
  x <- rt(N, df=20)
  x <- rexp(N)
  x <- c(rnorm(N/2), rnorm(N/2, sd=2))
#  x <- rnorm(N)

  c(shapiro.test(x)$p.value,
    ad.test(x)$p.value,
    cvm.test(x)$p.value,
    lillie.test(x)$p.value
  )
  })

rowMeans(res<0.05)
}

library("parallel")

mclapply(N, runsim, mc.cores=3L)
[[1]]
[1] 0.127 0.119 0.115 0.096

[[2]]
[1] 0.231 0.217 0.199 0.159

[[3]]
[1] 0.351 0.344 0.314 0.223

[[4]]
[1] 0.591 0.562 0.523 0.368

Wallyplot

asas

# simuler med var het
N <- 200
x <- rbinom(N, size=1, prob=.5)
y <- rnorm(N, mean=.5*x, sd=1+x)
t.test(y ~ x, var.equal=TRUE)

    Two Sample t-test

data:  y by x
t = -1.0557, df = 198, p-value = 0.2924
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.6176365  0.1869158
sample estimates:
mean in group 0 mean in group 1 
      0.1828423       0.3982026 
t.test(y ~ x, var.equal=FALSE)

    Welch Two Sample t-test

data:  y by x
t = -1.0297, df = 147.54, p-value = 0.3048
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.6286782  0.1979575
sample estimates:
mean in group 0 mean in group 1 
      0.1828423       0.3982026

Vis at den er robust

library("mlbench")
data(PimaIndiansDiabetes2) 

library("ggplot2")
data(midwest)
ggplot(midwest, aes(x=area, y=poptotal)) + geom_point()

library("nortest")
qqq <- sapply(1:5000, function(x) { 
  dd <- exp(rnorm(100, mean=5))
#  dd <- rt(100, df=20)
#  dd <- 
  c(shapiro.test(dd)$p.value,
    ad.test(dd)$p.value,
    cvm.test(dd)$p.value,
    lillie.test(dd)$p.value)
  } )

# pairs(t(qqq))
rowMeans(qqq<0.05)

Så starter det.

Tabel 1: aasdasdadsasdasd asd a asd ad
Girth Height Volume
8.3 70 10.3
8.6 65 10.3
8.8 63 10.2
10.5 72 16.4
10.7 81 18.8
10.8 83 19.7
11.0 66 15.6
11.0 75 18.2
11.1 80 22.6
11.2 75 19.9
11.3 79 24.2
11.4 76 21.0
11.4 76 21.4
11.7 69 21.3
12.0 75 19.1
12.9 74 22.2
12.9 85 33.8
13.3 86 27.4
13.7 71 25.7
13.8 64 24.9
14.0 78 34.5
14.2 80 31.7
14.5 74 36.3
16.0 72 38.3
16.3 77 42.6
17.3 81 55.4
17.5 82 55.7
17.9 80 58.3
18.0 80 51.5
18.0 80 51.0
20.6 87 77.0 
f <- function(n, veff=0, p=.25) {
  df <- data.frame(x=rbinom(n, size=1, prob=p)) %>% mutate(y=rnorm(n, mean=0, sd=(1+x*veff)))
  df
}

result <- sapply(1:5, function(ve) {
  res <- f(100, veff=ve) %>% lm(y ~ x, data=.) %>% summary() %>% coef() 
  res[2,2]
}
)
comments powered by Disqus