3 Deskriptiv statistik

3.1 Mean, median, SD and quantiles

Eksempel: Vitamin D (uge 1)

Se kapitel 6.2 i Descriptive analysis i R-intro (herunder brug af funktionen tapply til at lave beregninger gruppevist).

F.eks.

tapply( vit$vitd, vit$country, mean)

   DK    SF    EI    PL 
47.17 47.99 48.01 32.56 

Teknisk: Vil du have det i en fin tabel a la mine slides kan du benytte min funktion sumStat() (se evt kapitel om at definere egne funktioner). Kopier først denne funktion ind i R

sumStat <- function(x){
  ntotal=length(x)   
  nmiss=ntotal-table(is.na(x))[1] # antal manglende vaerdier
  n=ntotal-nmiss                  # antal med observerede vaerdier
  mean=mean(x,na.rm=T)
  median=median(x,na.rm=T)
  Q1=quantile(x, probs=.25,na.rm=T)
  Q3=quantile(x, probs=.75,na.rm=T)
  SD=sd(x,na.rm=T)
  min=min(x,na.rm=T)
  max=max(x,na.rm=T)
  samlet <- data.frame(n,nmiss,mean,median,SD, Q1,Q3,min,max)
  row.names(samlet) <- NULL
  return( samlet ) 
}

Dernæst kan vi udføre beregningerne enten for et land af gangen (her Irland)

sumStat( vitEI$vitd )

   n nmiss  mean median    SD   Q1   Q3 min   max
1 41     0 48.01   44.8 20.22 34.4 62.7  17 110.4

eller for alle 4 lande på en gang

tx  <- tapply(vit$vitd,vit$country,sumStat)
tx

$DK
   n nmiss  mean median    SD   Q1   Q3  min  max
1 53     0 47.17   47.8 22.78 27.6 61.8 11.4 93.6

$SF
   n nmiss  mean median    SD    Q1   Q3 min  max
1 54     0 47.99   46.6 18.72 35.52 60.9 5.2 96.6

$EI
   n nmiss  mean median    SD   Q1   Q3 min   max
1 41     0 48.01   44.8 20.22 34.4 62.7  17 110.4

$PL
   n nmiss  mean median    SD Q1   Q3 min  max
1 65     0 32.56   32.5 12.46 25 39.9 5.4 60.3

Teknisk: Vi kan klistre linjerne sammen ved at kalde rbind (RowBIND):

do.call('rbind', tx)

    n nmiss  mean median    SD    Q1   Q3  min   max
DK 53     0 47.17   47.8 22.78 27.60 61.8 11.4  93.6
SF 54     0 47.99   46.6 18.72 35.52 60.9  5.2  96.6
EI 41     0 48.01   44.8 20.22 34.40 62.7 17.0 110.4
PL 65     0 32.56   32.5 12.46 25.00 39.9  5.4  60.3

3.2 Fraktiler

Eksempel: Vitamin D (irske kvinder, uge 1)

2.5 og 97.5%-fraktiler

quantile( vitEI$vitd, probs = c(.025,.975) )

 2.5% 97.5% 
 18.0  89.1 

Vi kan evt tvinge R til at beregne dem på samme måde som i SAS/SPSS, hvor der benyttes en interpolation mellem de observerede værdier (ekstra argument type=6, der findes 9 forskellige typer at beregne fraktiler på i R …):

quantile( vitEI$vitd, probs = c(.025,.975), type=6 )

  2.5%  97.5% 
 17.05 109.33 

3.2.1 Bestemmelse af fraktiler i t-fordelingen

Hertil benyttes qt() (Quantiles T-distribution). Øvre 2.5%-fraktil med 40 frihedsgrader bestemmes som:

qt( 0.975, df=40 )

[1] 2.021

3.3 Tabellering

Eksempel: Vitamin D (irske kvinder, uge 1)

Ønsker vi at vide hvor mange der har e.g. Vitamin D niveau under 25 laver vi en tabel delt op på om Vitamin D er over eller under 25:

tab <- table( vitEI$vitd < 25 )
tab 

FALSE  TRUE 
   36     5 

Vi får en tabel med TRUE/FALSE, hvor TRUE angiver antallet med vitd<25 (her 5 kvinder).

Vi bestemmer andelen af kvinder med Vitamin D under 25 til 12%:

prop.table( tab )

FALSE  TRUE 
0.878 0.122 

3.4 Konfidensintervaller

3.4.1 CI for et gennemsnit

Beregnes enten i hånden ud fra formlen

\[ \bar Y \pm q_{0.975}\times \frac{\rm SD}{\sqrt{n}} \] hvor \(q_{0.975}\) er 97.5%-fraktilen i en t-fordeling med \(n-1\) frihedsgrader (se hvordan du bestemmer den i afsnit 3.2.1 ovenfor).

Eller den smarte løsning: Lav et one-sample t-test, så får du automatisk gennemsnit og CI ud (se kap 4.1)

3.4.2 CI for en spredning

Hertil benyttes pakken TeachingDemos som først installeres

# Denne linje skal kun køres en gang
install.packages("TeachingDemos")

og nu bestemmes CI for spredningen på vitamin D for de danske kvinder:

library(TeachingDemos)
sigma.test( vitDK$vitd )

    One sample Chi-squared test for variance

data:  vitDK$vitd
X-squared = 26991, df = 52, p-value <2e-16
alternative hypothesis: true variance is not equal to 1
95 percent confidence interval:
 365.7 794.6
sample estimates:
var of vitDK$vitd 
            519.1 

Vi får her variansen (spredningen\(^2\)) incl CI for, så vi skal tage kvadratrod for at få spredningen (se evt anvendelsen i forelæsningsscript uge 2).

Bemærk at vi yderst sjældent (som i aldrig!?) er interesseret i test af specifik værdi for spredningen (default er hypotesen \(H_0: \sigma=1,\) kan ændres med ekstra argument sigma=2 for at teste at spredningen er 2). Dvs i output ovenfor vil vi typisk kun aflæse CI.