Statistikpakken R
ST501 – R NOTER
Søren Møller og Søren Hesel
April 1, 2011
Dette er et kort notesæt, hvor der skulle være alt hvad der skal bruges i kurset ST501. Ønskes
der et mere udførligt notesæt henvises til http://cran.r-project.org/doc/manuals/R-intro.pdf
eller www.sbtc.ltd.uk/notes/Rintro.pdf. Det første er et generelt notesæt, mens det andet er
mere om vektorregning i R.
Vigtigt: R har en ganske god hjælpefunktion; ved man for eksempel ikke hvordan man får
længden af en vektor søger man ganske enkelt på længde (på engelsk selvfølgelig). Det vil sige
man skriver:
help(length)
For at komme ud af hjælp skrives q.
1
Installation af R
R kan downloades gratis fra hjemmesiden http://www.r-project.org til diverse styresystemer
(linket CRAN på venstre side), og installeres let. Det er installeret på maskinerne vi bruger til
computerøvelserne. I linux søges efter r-base i f.eks. Synaptic.
2
Om R
R er et vektor- og matrixbaseret program, hvilket vil sige alle variabler opfattes som en matrix.
Man tilskriver x en værdi med kommandoen <-. Dvs. et eksempel er
x <- 5
Vil man oprette en vektor f.eks. (1, 2, 3) skrives:
x <- c(1,2,3)
dvs. kommandoen c laver en vektor.
Skal man skrive til den i’te række og j’te søjle af en matrix x (f.eks. give den værdien 5)
skrives:
x[i,j] <- 5
og skal kun j’te søjle bruges:
1
x[,j] <- c(5,5,5)
Bemærk at R indekserer fra 1, hvilket vil sige, at når man tæller pladser i sine søjler starter
man fra 1.
Bemærk, dette også medfører, at når man f.eks. skal lave et plot af data eller en matrix skal
det specificeres om det kun skal være af en søjle eller vektor - ellers bruges hele matricen eller
datasættet som data.
2.1
Flere vektorfunktioner i R
• rep(a,b): Gentager a b gange: dvs. laver vektoreren ( ).
• Eksempel: rep(1,5) giver (1, 1, 1, 1, 1)
• seq(a,b,c): Laver en sekvens af tal fra a til b med skridtlængde c.
• Eksempel: seq(1,5,0.5) giver (10 15 20 25 30 35 40 45 50).
2.2
Data i R
R kan arbejde med næsten alle slags data-filer – ikke bare tekstfiler. Dog skal vi ikke bruge
andre slags i dette kursus. Det er vigtigt, at hvis man har overskrift på sin data (se nedenfor),
må der ikke være mellemrum i ordene. Decimaltal skal være med punktum som decimaldeler,
da R bruger angelsaksisk decimaltal ligesom Maple.
3
Opstart af R
I R skrives kommandoer i et R-script; dette er den bedste måde at kunne gengive det man har
lavet, da det ellers er svært at gemme, hvad man har lavet. En script-fil kan gemmes som speciel
R-fil, men en almindelig tekst-fil kan også bruges. I en scriptfil bruges # som udkommentering
– dvs. at R ikke læser noget efter # på en linje.
3.1
Linux
Under linux er R et konsol-program og kaldes med kommandoen R fra den mappe hvor ens
datafiler er. Man kører sine R-scripts med kommandoen source("filnavn"). For at udskrive tal
bruges kommandoen print.
3.2
Windows
R startes under Windows fra startmenuen. Det er vigtigt at man som første handling efter
opstart vælger Change dir i menuen File, og udvælger den mappe, som indeholder ens datafiler.
Man kører et script ved at åbne det med File, Open Script og så markere det, der skal køres og
Ctrl+R. Linux-måden kan også bruges.
2
4
4.1
Betjening af R
Linux
Kommandoer indtastes direkte i konsollen, plots fremkommer i et ekstra vindue.
4.2
Windows
Kommandoer indtastes i det vindue der hedder R Console, plots fremkommer i et ekstra vindue.
5
Regning med vektorer
Meget af R’s vektorregning er lige som med programmet Maple. Man lægger sammen og trækker
fra pladsvis med + og -. Man ganger og dividerer pladsvis med * og /, og skal man lave
vektorprodukt og matrixprodukt bruges %*%.
5.1
Transformation af variable
I statistik er der ofte en ikke-lineær sammenhæng mellem og så for eksempel sammenhængen
log ∼ ja så må man ud fra lave en vektor med værdierne af log . Dette gøres ved at
tilskrive en ny vektor værdierne log .
5.2
Eksempel
> x
[1] 1 2 3 3 4 5 6 8 9 11
> xx <- log(x)
> xx
[1] 0.0000000 0.6931472 1.0986123 1.0986123 1.3862944 1.6094379 1.7917595
[8] 2.0794415 2.1972246 2.3978953
Ofte skal man bruge tal for længden af en vektor, summen af tallene og så videre. Dette er
gjort meget intuitivt i R med henholdvis kommandoerne length(x) og sum(x).
6
Indlæsning af data
Data indlæses med kommandoen
data <- read.table(filename,header)
hvor filename er navnet på din datafil (som helst skal ligge i den mappe du valgte som aktiv
mappe, da du startede R) og header er header=FALSE når datafilen ikke indeholder variabelnavne hhv. header=TRUE når filen indeholder variabelnavne. Betegnelsen data er det objektnavn
vi giver til vores data, og som vi kan bruge for at arbejde med dem. I dette notesæt kalder vi
vores data for data, x eller y, men andre navne kan sagtens bruges.
De indlæste data kan udskrives på skærmen ved at indtaste variablenavnet som kommando.
3
data
Hvis datafilen indeholder mere end en søjle (f.eks. x- og y-værdier til en lineær regression),
kan disse separeres ved kommandofølgen
data <- read.table(filename,header)
x <- data$V1
y <- data$V2
Vær opmærksom på hvilken søjle der er hvilken, og udskiv om nødvendigt x og y i kommandofølgen.
6.1
Eksempel
> data <- read.table("data.txt",header=FALSE)
> data
V1
1 18.95
2 19.00
3 17.95
4 15.54
5 14.00
6 12.95
7
8.94
8
7.49
9
6.00
10 3.99
6.2
Eksempel
> data <- read.table("data2.txt",header=FALSE)
> data
V1
V2
1
1 18.95
2
2 19.00
3
3 17.95
4
3 15.54
5
4 14.00
6
5 12.95
7
6 8.94
8
8 7.49
9
9 6.00
10 11 3.99
> x<-data$V1
> y<-data$V2
> x
4
[1] 1 2 3 3 4 5 6 8 9 11
> y
[1] 18.95 19.00 17.95 15.54 14.00 12.95
7
8.94
7.49
6.00
3.99
Fordelinger
For at få værdier af tætheds- og fordelingsfunktioner bruges henholdsvis d og p foran den ønskede
fordelings navn. Vi laver et eksempel med den vigtigste fordeling, normalfordelingen, nedenfor.
Kvantilerne for en fordeling fås med bogstavet q.
Man skriver et af bogstaverne foran betegnelsen for den pågældende fordeling, variablen og til
sidst parametre til at finde den rigtige fordeling - se eksempel med normalfordelingen nedenfor.
7.1
Normalfordelingen
Normalfordelingen bruges tit som antagelse for fordelingen af data. Derfor er det vigtigt, at
kunne få tæthedsfunktion og fordelingsfunktion for denne fordeling. Efter det bogstav, der
repræsenterer hvilken funktion man ønsker skriver man norm. Derefter skal variablen skrives,
og til sidst middelværdi og spredning (sd), da disse to alene karakteriserer normalfordelingen.
Se afsnit 12 for eksempel for brug af kvantil-funktion.
7.1.1
Eksempel
> dnorm(0.5,0,1)
[1] 0.3520653
Fordelingsfunktionen fås ved pnorm.
7.1.2
Eksempel
> pnorm(0.5,0,1)
[1] 0.6914625
7.2
Andre fordelinger
De følgende udbredte fordelinger kan også bruges i R på tilsvarende måde:
Fordeling
Student -t
Fisher -F
2
Binomial
Kode
t
f
chisq
binom
Table 1: Flere fordelinger
Parametre
Antal frihedsgrader
Antal frihedsgrader k og l
Antal frihedsgrader
Populationsstørrelse og ssh. for succes.
5
8
Plots
8.1
Scatterplot
Skal der laves et plot af mod , bruges kommandoen plot.
8.1.1
Eksempel
> data
V1
V2
1 2.00 5.00
2 2.50 4.50
3 2.30 5.40
4 1.90 4.70
5 2.50 6.70
6 2.80 7.00
7 3.60 4.50
8 2.70 5.20
9 2.10 4.10
10 1.89 4.78
plot(data$V1,data$V2)
Figure 1: Scatterplot
8.2
Histogrammer
Histogrammer genereres med kommandoen
6
hist(data)
Som regel skal man lægge en fordeling oveni, og så kan man bede om et tæthedshistogram;
dette gøres ved
hist(data,probability=T)
Skal der indlægges en normalfordelingskurve oveni histogrammet, bruges i stedet kommandofølgen
m<-mean(data)
v<-var(data)
x<-seq(-5,5,.1)
n<-dnorm(x,m,sqrt(v))
lines(x,n)
hist(data)
Bemærk her, at det antages at vi skal bruge intervallet [−5 5]. Ellers skal der gengives andre
tal i tredje linie.
8.2.1
Eksempel
> data
[1] 1 2 3
> hist(data)
3
4
5
6
8
9 11
Figure 2: Histogram
7
8.3
Histogram med ikke-numeriske data
Det sker, at ens variabel for eksempel er antal gange en farve er observeret. Hvis man vil lave
et histogram for sådanne data bruges:
8.3.1
Eksempel
n <- length(data$farve)
plot(1:n,data$antal,type="h",xlab="",labels=F,tick=F)
axis(1,at=1:n,labels=data$farve,cex.axis=0.75,las=2)
axis(2,seq(1,floor(max(data$antal))+1,5),las=1)
8.4
Boxplots
Boxplots genereres med kommandoen
boxplot(data)
8.4.1
Eksempel
> data
[1] 1 2 3 3
> boxplot(data)
4
5
6
8
9 11
Figure 3: Boxplot
8
8.5
Kvantilplots
Normalfordelings-kvantilplots kan genereres med kommandofølgen
n<-length(data)
q<-qnorm((1:n)/(n+1))
plot(q,sort(data))
8.5.1
Eksempel
> data
[1] 1 2 3 3 4 5 6 8 9 11
> n<-length(data)
> n
[1] 10
> q<-qnorm((1:n)/(n+1))
> q
[1] -1.3351777 -0.9084579 -0.6045853 -0.3487557 -0.1141853
[7] 0.3487557 0.6045853 0.9084579 1.3351777
> plot(q,sort(data))
0.1141853
Kommandoen qqnorm kan bruges, når data er normalfordelt. Så skrives der
qqnorm(data)
Figure 4: Kvantilplot
9
9
Middelværdi og varians
Middelværdie af en datavektor bestemmes med kommandoen
mean(data)
og variansen bestemmes med
var(data)
Da spredningen er kvadratroden af variansen kan denne bestemmes ved
sqrt(var(data))
eller direkte ved kommandoen sd,
sd(data)
ønskes der et vægtet gennemsnit (f.eks. hvis antal målinger af en bestemt værdi er angivet
i en seperat søjle), bruges i stedet kommadoen
weighted.mean(data,weight)
9.1
Eksempler
> data
[1] 1 2 3 3 4
> mean(data)
[1] 5.2
> var(data)
[1] 10.62222
> sqrt(var(data))
[1] 3.259175
10
5
6
8
9 11
Median og kvartiler
Median, øvre og nedre kvartil samt mindste og største værdi i datasættet kan vises med kommandoen
summary(data)
10.1
Eksempel
> data <- read.table("data.txt",header=FALSE)
> data
V1
1 18.95
2 19.00
10
3 17.95
4 15.54
5 14.00
6 12.95
7
8.94
8
7.49
9
6.00
10 3.99
> summary(data)
V1
Min.
: 3.990
1st Qu.: 7.853
Median :13.475
Mean
:12.481
3rd Qu.:17.348
Max.
:19.000
Her angiver Min. mindste og Max. største måling, mens 1st Qu. og 3rd Qu. angiver nedre
og øvre kvartil, Median medianen og Mean middelværdien.
11
Lineær regression
For at lave lineær regression mellem to vektorer x og y bruges kommandoen
fit <- lm(y~x)
Så kan informationer om regressionen vises med kommandoerne
fit
summary(fit)
11.1
Eksempel
> data<-read.table("data2.txt",header=FALSE)
> data
V1
V2
1
1 18.95
2
2 19.00
3
3 17.95
4
3 15.54
5
4 14.00
6
5 12.95
7
6 8.94
8
8 7.49
9
9 6.00
11
10 11 3.99
> x<-data$V1
> y<-data$V2
> x
[1] 1 2 3 3 4 5 6 8 9 11
> y
[1] 18.95 19.00 17.95 15.54 14.00 12.95
> fit<-lm(y~x)
> fit
8.94
7.49
6.00
3.99
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept)
21.151
x
-1.667
Her er tallet under (Intercept) skæring med y-aksen og tallet under x hældningen for
bedste rette linje.
> summary(fit)
Call:
lm(formula = y ~ x)
Residuals:
Min
1Q Median
-2.2072 -0.5206 -0.2340
3Q
0.9183
Max
1.8010
Coefficients:
Estimate Std. Error t value
(Intercept) 21.1509
0.7558
27.98
x
-1.6673
0.1249 -13.35
--Signif. codes: 0 ’***’ 0.001 ’**’ 0.01
Pr(>|t|)
2.87e-09 ***
9.51e-07 ***
’*’ 0.05 ’.’ 0.1 ’ ’ 1
Residual standard error: 1.222 on 8 degrees of freedom
Multiple R-squared: 0.957,
Adjusted R-squared: 0.9516
F-statistic: 178.1 on 1 and 8 DF, p-value: 9.506e-07
Her ses skæring med y-aksen og hældning i søjlen under Estimate mens Multiple R-squared
angiver determinations-koefficienten 2 .
12
11.2
Anovatabeller
Når man har lavet en lineær model i R, gemt i objektet fit laves en Anovatabel med kommandoen anova(fit).
11.2.1
Eksempel
anova(fit)
Analysis of Variance Table
Response: data$V2
Df Sum Sq Mean Sq F value
Pr(>F)
data$V1
1 265.751 265.751 178.09 9.506e-07 ***
Residuals 8 11.938
1.492
--Signif. codes: 0 ´S***Š 0.001 ´S**Š 0.01 ´S*Š 0.05 ´S.Š 0.1 ´S Š 1
>
11.3
Konfidensinterval for lineær regression
For at lave et konfidensinterval for parametrene i en model bruges confint.
11.3.1
Eksempel
> confint(fit)
2.5 %
97.5 %
(Intercept) 19.407867 22.893848
x
-1.955388 -1.379173
12
Konfidensinterval for middelværdi
Hvis data antages at være normalfordelte med kendt varians kan man lave et niveau 1 −
konfidensinterval for middelværdien med formlen:
¡
¢
− SE()2 + SE()2
hvor er gennemsnittet af data, SE() = og 2 er 1 − 2-kvartilen for standard
normalfordelingen N(0 1), som i dette tilfælde fås med qnorm(1-alpha/2,0,1).
13
© Copyright 2024