I statistik vil vi ofte simulere data for at se, hvordan disse opfører sig under idealiserede forhold. Dette kan være noget så simpelt som at se, hvad der sker, når vi kaster en "fair" terning. Som et eksempel kan vi undersøge, hvor ofte summen af antal øjne på to terninger er større end eller lig med 9. Vi vil derfor gerne kunne få R til at kaste terningerne for os. Man taler om at generere tilfældige tal (random number generator). Følgende tabel viser nogle af mulighederne i R. $$\begin{array}{ll} \hline \text{Navn} & \text{Beskrivelse} \\ \hline \text{sample(6,10,replace=TRUE)} & \text{10 tilfældige hele tal mellem 1 og 6} \\ \text{rbinom(10,20,0.5)} & \text{10 tilfældige udfald fra en binom(20,0.5)} \\ \text{rpois(10,3)} & \text{10 tilfældige udfald fra en pois(3)} \\ \text{rmultinom(10,3,c(0.1,0.6,0.3))} & \text{10 tilfældige udfald fra en} \\ & \text{multinom(3,(0.1,0.6,0.3))} \\ \text{rnorm(10,mu,sigma)} & \text{10 tilfældige udfald fra en N(mu,sigma)} \\ \text{randperm(6)} & \text{tilfældig permutation af tallene 1 til 6} \\ \hline \end{array}$$ Som et første eksempel vil jeg svare på det spørgsmål, der blev stillet ovenfor i forbindelse med kast af to terninger. Jeg vil kaste de to terninger 10$\,$000 gange og se, hvor stor en andel af disse hvor summen af antal øjne er større end eller lig med 9. De 10$\,$000 kast af den ene terning genereres med sample(6,10000,replace=T), hvor "6" angiver, at de mulige udfald er $1,2,3,4,5,6$ svarende til en terning, og "10000" siger, at der skal laves en vektor med 10$\,$000 tilfældige udfald.

 > x=sample(6,10000,replace=T) # 1.terning kastes 10000 gange
 > y=sample(6,10000,replace=T) # 2.terning kastes 10000 gange
 > z=x+y                       # sum af antal øjne i de 10000 kast
 > sum(z>=9)/10000             # andel med sum >=9

Som et andet eksempel vil jeg simulere $t$-teststørrelsen fra afsnit 4.4 og sammenligne med $t$-fordelingen. Teststørrelsen er på formen $t=\bar x/\sqrt{s^2/n},$ hvor $n$ er antallet af observationer, $\bar x$ er gennemsnittet af disse, og $s^2$ er den empiriske varians.

 > n=10
 > nSim=1000000
 > xbar=rnorm(nSim,0,1/sqrt(n))
 > s2=rchisq(nSim,n-1)/(n-1)
 > t=xbar/sqrt(s2)

 > sum(abs(t)>qt(0.975,n-1))/Nsim
 [1] 0.050297

I R kan man finde fordelingsfunktionen, det vil sige sandsynligheden for at ligge til venstre for et punkt, for alle standard fordelinger. Disse er navngivet ved bogstavet p (p står for probability) efterfulgt af navnet på fordelingen. De vigtigste kan ses i listen herefter. De tilsvarende funktioner for at finde fraktiler i en fordeling fremkommer ved at erstatte det indledende p med q (q står for quantile). Ved en standard normalfordeling kan man udelade mu=0 og sigma=1 i kaldet til pnorm. $$\begin{array}{ll} \hline \text{Navn} & \text{Beskrivelse} \\ \hline \text{pnorm(1.96,mu,sigma)} & N_{\text{cdf}}(1.96,\text{mu},\text{sigma}) \\ \text{pbinom(3,10,0.5)} & P(X\leq 3)\enspace X\sim\text{binom(10,0.5)} \\ \text{ppois(3,1.5)} & P(X\leq 3),\enspace X\sim\text{pois(1.5)} \\ \text{pchisq(5.99,2)} & \chi_{\text{cdf}}(5.99,2) \\ \text{pt(4.3,2)} & t_{\text{cdf}}(4.3,2) \\ \text{pf(2.9,7,12)} & F_{\text{cdf}}(2.9,7,12) \\ \hline \end{array}$$

 > pnorm(1.96,0,1)
  [1] 0.9750021

 > pnorm(1.96)
  [1] 0.9750021

 > pchisq(5.99,2)
  [1] 0.9499634

 > qchisq(0.95,2)
  [1] 5.9914645

For en vektor $x$ med observerede dataværdier beregner man ofte forskellige beskrivende mål. Tabellen nedenfor gengiver de mest almindelige. $$\begin{array}{ll} \hline \text{Navn} & \text{Beskrivelse} \\ \hline \text{mean(x)} & \text{gennemsnit} \\ \text{var(x)} & \text{empirisk varians} \\ \text{sd(x)} & \text{empirisk spredning} \\ \text{median(x)} & \text{median} \\ \hline \end{array}$$ For at få en direkte fornemmelse af fordelingen kan man lave et histogram. I R gøres dette med funktionen hist. Denne funktion laver som default en histogramfigur og danner en struktur, hvorfra antallene i de forskellige intervaller kan findes. Hvis $x$ er en vektor med dataværdierne, kan man enten skrive

hist(x)$\quad$eller$\quad$hist(x,breaks=Nint)$\quad$eller$\quad$hist(x,breaks=endePkt)

I den første version beslutter R sig selv for en inddeling af aksen, der dækker alle observationerne. I den anden version specificerer man antallet Nint af intervaller, og i den sidste version angiver man en vektor endePkt med intervalendepunkter (alle værdierne i $x$ skal være indeholdt i de angivne intervaller, ellers får man en fejlmeddelelse). For et interval fra $a$ til $b$ tæller R antallet af observationer $x_i$ med $a< x_i\leq b.$ I de to første versioner, samt i den tredje hvis intervallerne defineret gennem endePkt alle har samme længde, laver R et antalshistogram, hvor højden på kasserne angiver antallet af observationer i intervallerne. For at få et tæthedshistogram, således at det samlede areal under kasserne er 1, skal man tilføje probability=TRUE til kaldet af hist. Man kan sætte sine egne titler på akserne og figuren ved brug af xlab, ylab og main som i beskrivelsen af plot. For at få antallene i hvert interval kan man skrive

antal=hist(x,breaks=endePkt)$\text{\textdollar}$counts

hvorefter antallene ligger i vektoren antal. Et jordstykke er delt op i 100 lige store områder. I hvert af disse områder er der plantet 50 ærter, og der registreres hvor mange af disse, der spirer og vokser op. De målte antal betegnes $x_i,$ $i=1,\ldots,100.$ Nedenfor laves et histogram med intervallængde 2 og efterfølgende et goodness of fit test for, om data kan beskrives med en binomialfordeling. Metoden er omtalt i afsnit 3.4

 > x =c(31,34,31,35,37,35,37,25,32,35,33,40,31,34,30,37,41,34,36,
   31,28,33,40,28,32,31,39,34,38,34,35,36,36,32,36,35,38,40,
   37,42,38,37,33,34,35,36,37,32,36,39,32,37,28,34,33,31,36,
   41,38,32,33,39,34,37,36,40,32,39,39,38,34,34,38,35,31,33,
   32,28,38,37,34,36,35,37,34,37,34,40,39,38,34,28,30,29,32,
   39,36,33,36,31)
 
 > n=length(x)  # antal observationer
 > phat=sum(x)/(50*n)  # skøn over p
 > endePkt=24+c(0:10)*2-0.5 # intervalendepunkter i histogram
 > antal=hist(x,breaks=endePkt)$counts  # figur med antalshistogram

 # beregne forventede fra binomialfordeling
 > p0=pbinom(endePkt,50,phat) # sandsynlighed til venstre for endepunkter
 > p=c(p0[2],p0[3:10]-p0[2:9],1-p0[10])  # intervalsandsynligheder
 # første interval fra 0 til og med 25, sidste fra 42 til 50

 > ex=n*p  # forventede
 > cbind(antal,round(ex,2)) # round afrunder tallene
      antal      
  [1,]     1  0.29
  [2,]     0  1.15
  [3,]     6  3.97
  [4,]    10 10.03
  [5,]    16 18.32
  [6,]    22 23.87
  [7,]    22 21.77
  [8,]    15 13.52
  [9,]     7  5.49
 [10,]     1  1.59
  
 # slår intervaller sammen for at få forventede >= 5
 > antal1=c(antal[1]+antal[2]+antal[3],antal[4:8],antal[9]+antal[10])
 > ex1=c(ex[1]+ex[2]+ex[3],ex[4:8],ex[9]+ex[10])
 > cbind(antal1,round(ex1,2))
     antal1      
 [1,]      7  5.41
 [2,]     10 10.03
 [3,]     16 18.32
 [4,]     22 23.87
 [5,]     22 21.77
 [6,]     15 13.52
 [7,]      8  7.08

 > gTest=2*sum(antal1*log(antal1/ex1));
 > df=length(antal1)-1-1;
 > pval=1-pchisq(gTest,df);

 > c(gTest,df,pval)
  [1] 1.1580229 5.0000000 0.9488405

I dette afsnit nævnes nogle af de funktioner i R, der bruges til klassiske statistiske tests. Hvis $x$ er en vektor med dataværdier, kan man lave et $t$-test for, at middelværdien har en bestemt værdi mu0 med t.test. Dette er gennemgået i afsnit 4.4. Den mest simple version er som følger:

t.test(x,mu=mu0)

Hvis mu=mu0 udelades, laves et test for, om middelværdien er nul. Hvis output placeres i testUD kan man adressere de forskellige dele af output ved at skrive testUD$\text{\textdollar}$navn. De vigtigste muligheder er: $$\begin{array}{ll} \hline \text{Navn} & \text{Beskrivelse} \\ \hline \text{statistic} & \text{værdien af t-teststørrelsen} \\ \text{parameter} & \text{antallet af frihedsgrader} \\ \text{p.value} & p\text{-værdien} \\ \text{conf.int} & \text{95\%-konfidensinterval for middelværdien} \\ \hline \end{array}$$ Konfidensintervallet er for den ukendte middelværdi, og størrelsen af konfidensintervallet bestemmes af et niveau, der kan vælges ved at tilføje conf.level=0.95 i kaldet til t.test. Funktionen t.test laver et såkaldt tosidet test, det vil sige, at både store positive og store negative værdier er kritiske. Hvis man ønsker et test, hvor det kun er store positive værdier, der er kritiske, skal man tilføje alternative="greater" i kaldet til t.test.Ovenfor blev $t$-testet for en observationsrække omtalt. Lad os nu betragte $t$-testet for samme middelværdi i to observationsrækker. Dette svarer til afsnit 4.13. Antag, at vektoren $x$ indeholder dataværdierne for den ene observationsrække og $y$ dataværdierne for den anden observationsrække. Den relevante kommando er

t.test(x,y,var.equal=TRUE)

som laver $t$-testet under antagelse om, at varianserne i de to grupper er ens. Hvis varianserne ikke er ens, skal man ændre TRUE til FALSE. Outputtet fra funktionen indeholder de samme elementer som ovenfor. Konfidensintervallet er i dette tilfælde for forskel i middelværdi mellem de to grupper $\delta=\mu_1-\mu_2.$I tilfældet med to grupper af observationer kan man lave et test for, at de to varianser er ens ved

var.test(x,y)

Igen vil outputtet have de komponenter, der er nævnt ovenfor. Konfidensintervallet er for forholdet mellem de to varianser. Output fra en funktion som t.test kaldes i R for en liste. Man kan se, hvilke elementer listen indeholder ved at skrive

names(navn)

hvor navn er navnet på den struktur output er placeret i. Den lineære regressionsmodel fra afsnit 5.1 analyseres i R ved hjælp af kommandoen lm(x$\sim$t). I dette udtryk er $x$ responsvariabel, og $t$ er den forklarende variabel. Udtrykket "x$\sim$t" kaldes modelformlen. Den følgende kode viser dette med Hubbles data omtalt i afsnit 5.6. De to figurer der produceres, kan ses i den sidste figur i afsnit R.2 øverste højre hjørne og nederste delfigur. En meget brugbar efterbehandling af resultatet fra lm er funktionen summary. Denne giver blandt andet skønnet over spredningen på data, og den giver en tabel med parameterestimater, standard error for disse og et $t$-test for, om parameteren kan antages at være nul. En anden efterbehandling er ordren confint(), der giver konfidensintervaller for parametrene i middelværdimodellen.

 > a=c(0.032,0.034,0.214,0.263,0.275,0.275,0.450,0.500,
    0.500,0.630,0.800,0.900,0.900,0.900,0.900,1.000,
    1.100,1.100,1.400,1.700,2.000,2.000,2.000,2.000)
 > h=c(170,290,-130,-70,-185,-220,200,290,270,200,300,
    -30,650,150,500,920,450,500,500,960,500,850,800,1090)
 > plot(a,h,xlab="Afstand",ylab="Hastighed", main="Regression")

 > lmUD=lm(h~a)  # model analyseres og resultat placeres i lmUD

 > lmUD$coefficients  # parameterestimater
  (Intercept)           a 
    -40.78365   454.15844 

 > confint(lmUD)  # konfidensintervaller
                  2.5 %   97.5 %
  (Intercept) -213.8253 132.2580
  a            298.1262 610.1906 

 > abline(lmUD)  # linje indtegnes

 > sumUD=summary(lmUD)  # uddrager information om den fittede model

 > sumUD
  Call:
  lm(formula = h ~ a)

  Residuals:
      Min      1Q  Median      3Q     Max 
  -397.96 -158.10  -13.16  148.09  506.63 

  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept)   -40.78      83.44  -0.489     0.63    
  a             454.16      75.24   6.036 4.48e-06 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Residual standard error: 232.9 on 22 degrees of freedom
  Multiple R-squared:  0.6235,	Adjusted R-squared:  0.6064 
  F-statistic: 36.44 on 1 and 22 DF,  p-value: 4.477e-06

 > df=sumUD$df[2]
 > s=sumUD$sigma
 > ssd=df*s^2
 > c(df,ssd,s,s^2)
  [1]      22.0000 1193442.3663     232.9107   54247.3803

 > plot(a,sumUD$residuals,xlab="Afstand",ylab="Residualer")
 > abline(0,0)

Funktionen lm kan også bruges til den generelle lineære model beskrevet i kapitel 6 og kapitel 7. Følgende eksempel ser på data i afsnit 6.6. For at beskrive modellen skal nogle af variablene defineres som faktorer.

 > areal=c(264,200,225,268,215,241,232,256,229,288,253,288,230,
   235,188,195,205,212,214,182,215,272,163,230,255,202,
   314,320,310,340,299,268,345,271,285,309,337,282,273,
   283,312,291,259,216,201,267,326,241,291,269,282,257)
 > stress=factor(rep(c("Uden","Med","Uden","Med"),c(13,13,13,13)))
 > lys=factor(rep(c("Lav","Hoej"),c(26,26)))

 > lmUD=lm(areal~lys*stress)  # tosidet variansanalyse
 
 > summary(lmUD)  # information om fittet til data

 Call:
 lm(formula = areal ~ lys * stress)

 Residuals:
     Min      1Q  Median      3Q     Max 
 -67.846 -19.385  -0.385  19.538  59.077 

 Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
 (Intercept)        268.846      8.299  32.395  < 2e-16 ***
 lysLav             -55.923     11.736  -4.765 1.79e-05 ***
 stressUden          35.231     11.736   3.002  0.00425 ** 
 lysLav:stressUden   -2.846     16.598  -0.171  0.86457    
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Residual standard error: 29.92 on 48 degrees of freedom
 Multiple R-squared:  0.5729,	Adjusted R-squared:  0.5462 
 F-statistic: 21.46 on 3 and 48 DF,  p-value: 5.867e-09
 
 # Lave F-test for additiv struktur
 > anova(lm(areal~lys+stress),lm(areal~lys*stress))

 Analysis of Variance Table

 Model 1: areal ~ lys + stress
 Model 2: areal ~ lys * stress
   Res.Df   RSS Df Sum of Sq      F Pr(>F)
 1     49 43003                           
 2     48 42976  1    26.327 0.0294 0.8646

Foregående