Erste Schritte mit R für Praktiker in den Biowissenschaften
V. Kiefel
04. März 2018
Contents
1 Hinweise zum Nachvollziehen der Beispiele
2 Ein-Stichproben-Analyse
2.1 Perzentilen, Histogramm
2.2 Mittelwert, Median und Standardabweichung
3 Einfache lineare Regression, Korrelationskoeffizient, Daten von einer Datei einlesen
4 Messung der Übereinstimmung der Beurteilung
4.1 Cohens Kappa
4.2 Spearmans Rangkorrelationskoefizient
5 Kontingenztafeln
5.1 Chiqudadrat-Test
5.2 Prüfung von k×2-Feldertafeln auf Trend
5.3 `Fisher's exakter Test '
5.4 Odds ratio
6 Vertrauensgrenzen relativer Häufigkeiten bei binominalverteilter Grundgesamtheit
7 Vergleich mehrerer Gruppen
7.1 Einfache Varianzanalyse
7.2 Bartlett Test
8 `Goodness of fit '-Tests (Anpassungstests an bekannte Verteilungen)
8.1 Kolmogorov-Smirnov Test für die Güte der Anpassung
8.2 Shapiro-Wilk-Test
8.3 Überprüfung von Genotyp-Häufigkeiten auf Vorliegen eines Hardy-Weinberg-Gleichgewichts
9 Vergleich von zwei unabhängigen Stichproben mit dem t-Test
10 Vergleich zweier abhängiger Stichproben im t-Test
11 Wilcoxon-Test für unabhängige Stichproben
12 Friedman-Test
13 H-Test (Kruskal-Wallis)
14 Überlebenszeit-Analyse
14.1 Beispieldaten aus der R-Implementation benutzen
14.2 Kaplan Meier-Kurve, Daten aus eigener Datentabelle einlesen
15 Statistische Funktionen
15.1 Chiquadrat-Verteilung
15.2 F-Verteilung
15.3 t-Verteilung (Student)
15.4 Standardnormalverteilung
16 Bestimmung einer nowendigen Stichprobengrö ße
16.1 Vergleich zweier Mittelwerte
16.2 Vergleich zweier relativer Häufigkeiten
17 Metaanalyse
18 Erstellung von Grafiken
18.1 Säulendiagramme
18.2 Farben und andere Grafikparameter
18.3 Boxplots
18.4 Scatterplots
18.5 Stripchart (eindimensionaler Scatterplot)
18.6 Linien verbinden Werte von beliebig vielen Individuen, die zwei oder mehr Zeitpunkten oder Wiederholungen entsprechen
18.6.1 Allgemeine Konstruktion mit plot()
18.6.2 Konstruktion ähnlicher Diagramme mit interaction.plot()
19 Verschiedenes
19.1 Berechnung der Determinante einer Matrix
19.2 Datumsberechnungen
20 Anhang
20.1 Hinweise für R unter Linux
20.1.1 Installation
20.1.2 R konfigurieren
20.1.3 PDF-, PostScript- und EPS-Dateien von Abbildungen erstellen.
20.2 Allgemeine Befehle
20.3 Einstellen von Optionen des Systems
20.4 Textdarstellung von R-Objekten als Textdatei speichern und wieder einlesen
20.5 Datendateien für R formatieren
20.6 Quellcodes
20.6.1 ToR.awk
20.6.2 Perc.R
21 Versionen
Index
Vorwort
Das vorliegende Manuskript ist eine kleine willkürlich zusammengestellte
Sammlung von ,,Kochrezepten'' für die Anwendung von R1
mit Verfahren, wie sie von Biowissenschaftlern und Medizinern ständig
benötigt werden.
Dieser Text wurde in der Absicht geschrieben, dem Benutzer die allerersten
Schritte mit R anhand von ihm vertrauten Problemen zu erleichtern.
Vor allem benötigt man häufig in der Eile ein Beispiel, anhand dessen man
rasch die Eingabe der Daten und die Formulierung einer Auswertung ausführen
kann. In diesem Sinne dient dieses Manuskript auch dem Autor als
Gedächtnisstütze.
Es empfiehlt sich dringend, parallel zum Ausprobieren der
Beispiele zu wichtigen
Funktionen die R-Dokumentation zu lesen, um zu erkennen, wie
die in diesem Text verwendeten Lösungen in der oft sehr knappen
Beschreibung ,,aussehen''. So wird man dann lernen, alternative und
möglicherweise bessere Lösungen zu finden.
Das Manuskript ersetzt nicht die ausgezeichnete Dokumentation zur
Beschreibung von R und es ersetzt nicht Lehrbücher oder
Originalarbeiten zu statischen Methoden. Vor oder gleichzeitig mit diesem
Dokument sollte man sich mit elementaren Datenstrukturen und Funktionen
vertraut machen, darunter Vektoren, Listen, Arrays, ,,data frames'', read.table(), c(), z. B. in dem Dokument ,,An introduction to
R'' des R development core teams, das jeder
Installation von R beigefügt ist [1].
Installation und Nutzung der Windows-Version von R dürften weitgehend
selbsterklärend sein. Hinweise auf Besonderheiten von R unter Linux
finden sich in Abschnitt 20.1.
Für die Richtigkeit der in diesem kleinen Manuskript gemachten Angaben wird
keine Gewähr übernommen. In fast allen Fällen ist der beschriebene Weg zur
Dateneingabe nicht der einzig mögliche. Bei zunehmender Vertrautheit mit
R sollte die Originaldokumentation und Fachliteratur zu den
zugrundeliegenden statistischen Verfahren zu Rate gezogen werden!
Empfehlenswerte Texte für den Einstieg in elementare Verfahren anhand von
R sind [1,2,3].
Kritzmow bei Rostock, im März 2018
V.K.
1 Hinweise zum Nachvollziehen der Beispiele
Die Nach Starten der R-Konsole können die Eingaben aus
dem Manuskript eingetippt werden, sie können aber auch aus dem HTML-Dokument
herauskopiert werden und an der Eingabeaufforderung der R Konsole mit
einem paste-Befehl eingefügt werden.
Kurze Eingaben sind im Folgenden durch die mit abgebildete Eingabeaufforderung:
>
am Zeilenbeginn gekennzeichnet. Bei langen Eingabezeilen wird durch das
+ in der Folgezeile die Fortsetzung einer langen Zeile signalisiert:
>
+
Bei der Eingabe sollen diese Zeichen > und + bei einer
zeilenweisen Eingabe nicht mit eingegeben oder hineinkopiert werden. Nach der
Eingabe, die mit Enter bestätigt wird, folgt dann die Ausgabe von
R ohne Eingabeaufforderung am Zeilenbeginn, nach der letzten Ausgabe
erscheint dann wieder die Eingabeaufforderung.
Eine alternative Form der Eingabe von R-Befehlen besteht darin,
R Code in eine Datei
zu schreiben und diese dann mit dem Befehl
> source("rcodedatei.R")
einzulesen. Ein Beispiel hierzu findet sich in dem Abschnitt
2.1. Hier werden Daten mit der Funktion c() in
einen Vektor eingefügt. Umfangreichere Datensätze schreibt man besser in eine
Tabelle (zum Beispiel mit einem Texteditor) und liest sie von R mit
der Funktion read.table() ein, ein Beispiel findet sich in Abschnitt
3.
2 Ein-Stichproben-Analyse
2.1 Perzentilen, Histogramm
Zur Berechnung von Perzentilen in der von Reed beschriebenen Form
[4] verwende man die Funktion perc in Perc.R
(Quellcode in Abschnitt 20.6.1). Dieser Abschnitt demonstriert die
Entwicklung und Verwendung eigener Funktionen.
# Aus Reed
# Anfang Eingabe
hapto <-
c(14, 21, 21, 30, 32, 35, 36, 36, 40, 44, 47, 48, 48, 48, 50, 51, 52, 54, 58,
59, 59, 62, 65, 66, 67, 67, 69, 71, 72, 76, 76, 77, 77, 77, 77, 78, 79, 79,
80, 81, 82, 84, 85, 86, 87, 87, 88, 88, 89, 90, 90, 93, 94, 94, 95, 96, 96,
97, 98, 98, 100, 100, 101, 101, 101, 103, 105, 106, 108, 108, 108, 109, 113,
114, 114, 114, 116, 116, 119, 126, 128, 129, 135, 136, 141, 142, 147, 147,
150, 161, 162, 170, 174, 174, 176, 179, 181, 191, 199, 225)
# Ende Eingabe
Wenn die Zuweisung der Werte für hapto in der Datei
hapto.R gespeichert ist und das Skript mit der Funktion
perc() (Abschnitt 20.6.2)
in der Datei perc.R, dann können Daten und Funktion
mit der Funktion source() geladen werden und die benutzerdefinierte
Funktion perc() aufgerufen werden, gefragt sei nach der
97.5-Perzentile:
> source("hapto.R")
> source("perc.R")
> perc(hapto,97.5)
[1] 194.8
>
Die eingebaute Funktion von R ergibt einen etwas anderen Wert:
> quantile(hapto, 0.975)
97.5%
186.25
>
wenn mehr als ein Percentil (Quantil) abgefragt werden soll:
> quantile(hapto, c(0.025,0.5,0.975))
2.5% 50% 97.5%
25.275 90.000 186.250
>
Ein Histogramm kann mit dem Kommando gezeichnet werden:
hist(hapto)
Mit eigenen Achsenbeschriftungen, Achsengrenzen versehen:
> hist(hapto, xlim=c(0,250), ylim=c(0,25),
+ main="Distribution of Haptoglobin Concentration",
+ xlab="Haptoglobin Concentration", ylab="Frequency")
>
2.2 Mittelwert, Median und Standardabweichung
Mittelwert, Median und Standardabweichung sollen zu den Werten in
perc berechnet werden:
> source("hapto.R")
> mean(hapto)
[1] 95.25
> sd(hapto)
[1] 42.71786
> median(hapto)
[1] 90
3 Einfache lineare Regression, Korrelationskoeffizient,
Daten von einer Datei einlesen
Es soll die Ausgleichsgerade zu folgenden x- und y- Werten
berechnet werden
(das Beispiel stammt aus [5,Seite 103]):
XWERT YWERT
1 13 12
2 17 17
3 10 11
4 17 13
5 20 16
6 11 14
7 15 15
Die Daten werden in dieser Form in eine Textdatei mit dem Namen 1.dat
mit einem Texteditor geschrieben.
Sie werden anschliessend an der R-Eingabeaufforderung
so eingelesen und ausgewertet:
> tafel <- read.table("1.dat", header=TRUE)
> x <- tafel$XWERT
> y <- tafel$YWERT
> res <- lsfit(x, y)
> ls.print(res)
Residual Standard Error=1.6695
R-Square=0.5023
F-statistic (df=1, 5)=5.0463
p-value=0.0746
Estimate Std.Err t-value Pr(>|t|)
Intercept 7.7288 2.8621 2.7004 0.0428
X 0.4262 0.1897 2.2464 0.0746
Die nach diesem Verfahren geschätzte Ausgleichsgerade ist damit: y = 7,7288+ 0,4262.
Zu diesem Ergebnis kommt man auch mit der angemesseneren Funktion lm():
> lm(y ~ x)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
7.7288 0.4262
Der Betrag des Korrelationskoeffizienten kann als Quadratwurzel aus
R-square bestimmt werden oder direkt mit der Funktion
cor.test():
> cor.test(x,y)
Pearson's product-moment correlation
data: x and y
t = 2.2464, df = 5, p-value = 0.07461
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.09505521 0.95310396
sample estimates:
cor
0.7087357
4 Messung der Übereinstimmung der Beurteilung
4.1 Cohens Kappa
Daten aus
[6,Seite 458]:
> library(vcd)
> daten <- c(53,5,2,11,14,5,1,6,3)
> mdaten <- matrix(daten, nrow=3, byrow=T)
> mdaten
[,1] [,2] [,3]
[1,] 53 5 2
[2,] 11 14 5
[3,] 1 6 3
> Kappa(mdaten)
value ASE lwr upr
Unweighted 0.4285714 0.08728716 0.2574917 0.5996511
Weighted 0.4923077 0.10057994 0.2951746 0.6894407
>
Resultat in [6]: 0,429
4.2 Spearmans Rangkorrelationskoefizient
Beispiel aus [7,Seite 310]:
> a <- c(7,6,3,8,2,10,4,1,5,9)
> b <- c(8,4,5,9,1,7,3,2,6,10)
> cor.test(a,b, method="spearman")
Spearman's rank correlation rho
data: a and b
S = 24, p-value = 0.003077
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.8545455
5 Kontingenztafeln
5.1 Chiqudadrat-Test
Einzugeben sei die Tafel:
14 22 32
18 16 8
8 2 0
Die Werte in einen Vektor und anschließend in eine Matrix überführen:
> r <- c(14, 18, 8, 22, 16, 2, 32, 8, 0)
> rm <- matrix(r, nrow=3)
> rm
[,1] [,2] [,3]
[1,] 14 22 32
[2,] 18 16 8
[3,] 8 2 0
Berechnung mit der Funktion chisq.test
> chisq.test(rm)
Pearson's Chi-squared test
data: rm
X-squared = 21.5765, df = 4, p-value = 0.0002433
Warning message:
Chi-squared approximation may be incorrect in: chisq.test(rm)
>
Einzugeben sei die Vierfeldertafel [7,S. 270]:
14 85
4 77
> r <- c(15,85,4,77)
> mr <- matrix(r,byrow=T,nrow=2)
> mr
[,1] [,2]
[1,] 15 85
[2,] 4 77
> chisq.test(mr)
Pearson's Chi-squared test with Yates' continuity correction
data: mr
X-squared = 3.8107, df = 1, p-value = 0.05093
> chisq.test(mr,correct=F)
Pearson's Chi-squared test
data: mr
X-squared = 4.8221, df = 1, p-value = 0.02810
>
Bei Vierfeldertafeln ist offenbar als Vorgabe Yates' Kontinuitätskorrektur
eingestellt, was im zweiten Aufruf abgewählt wurde.
5.2 Prüfung von k×2-Feldertafeln auf Trend
Die folgende Tafel soll auf Trend untersucht werden
[7,Seite 365]:
geheilt
Therapie nach
Erfolg spez. Anzahl
Score Therapie therapiert
1 2 10
2 16 34
3 22 36
Eingabe [2]:
> geheilt <- c(2,16,22)
> summe <- c(10,34,36)
> prop.trend.test(geheilt,summe,score=c(1,2,3))
Chi-squared Test for Trend in Proportions
data: geheilt out of summe ,
using scores: 1 2 3
X-squared = 5.2197, df = 1, p-value = 0.02233
Die gleiche Tafel würde beim Vergleich ohne Berücksichtigungs des Trends
keine signifikanten Unterschiede erkennen lassen:
> tafel <- c(14, 18, 8, 22, 16, 2)
> tmatrix <- matrix(tafel, byrow = T, nrow = 2)
> tmatrix
[,1] [,2] [,3]
[1,] 14 18 8
[2,] 22 16 2
> chisq.test(tmatrix)
Pearson's Chi-squared test
data: tmatrix
X-squared = 5.4954, df = 2, p-value = 0.06407
5.3 ,,Fisher's exakter Test''
Einzugeben sei die Tafel:
10 4
2 8
> tf <- c(10, 2, 4 ,8)
> mtf <- matrix(tf, nrow=2)
> fisher.test(mtf, alternative="greater")
Fisher's Exact Test for Count Data
data: mtf
p-value = 0.01804
alternative hypothesis: true odds ratio is greater than 1
95 percent confidence interval:
1.435748 Inf
sample estimates:
odds ratio
8.913675
Hier wurde der p-Wert für eine einseitige Fragestellung ermittelt. Für
die zweiseitige Fragestellung:
> fisher.test(mtf, alternative="two.sided")
Fisher's Exact Test for Count Data
data: mtf
p-value = 0.03607
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.114920 123.433541
sample estimates:
odds ratio
8.913675
5.4 Odds ratio
Eine weitere Möglichkeit zur Berechnung von Konfidenzintervallen einer Odds
Ratio besteht in der Verwendung der Funktion oddsratio im Paket
epitools (Beispiel aus [8,Seite 490])
> library(epitools)
> tabelle <- matrix(c(24,96,48,592),nrow=2,byrow=T)
> tabelle
> oddsratio(tabelle,method="midp",conf.level=0.95)
$data
Outcome
Predictor Disease1 Disease2 Total
Exposed1 24 96 120
Exposed2 48 592 640
Total 72 688 760
$measure
odds ratio with 95% C.I.
Predictor estimate lower upper
Exposed1 1.000000 NA NA
Exposed2 3.085417 1.779715 5.236565
$p.value
two-sided
Predictor midp.exact fisher.exact chi.square
Exposed1 NA NA NA
Exposed2 0.0001004225 0.0001119003 1.780411e-05
$correction
[1] FALSE
attr(,"method")
[1] "median-unbiased estimate & mid-p exact CI"
>
6 Vertrauensgrenzen relativer Häufigkeiten bei
binominalverteilter Grundgesamtheit
Fragestellung: Man bestimme den 95%-Vertrauensbereich für 7
Treffer unter 20 Beobachtungen.
> binom.test(7,20, conf.level=0.95)
Exact binomial test
data: 7 and 20
number of successes = 7, number of trials = 20, p-value = 0.2632
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.1539092 0.5921885
sample estimates:
probability of success
0.35
>
Eine Alternative zur Eingabe wäre:
> binom.test(c(7, 13), conf.level=0.95)
Gelegentlich können große Zahlen nicht berechnet werden, dann kann als
Alternative prop.test verwendet werden.
> prop.test(7,20, conf.level=0.95)
1-sample proportions test with continuity correction
data: 7 out of 20, null probability 0.5
X-squared = 1.25, df = 1, p-value = 0.2636
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
0.1630867 0.5905104
sample estimates:
p
0.35
Als weitere Alternative steht die Funktion binom.exact() aus dem
Paket epitools zur Verfügung.
7 Vergleich mehrerer Gruppen
7.1 Einfache Varianzanalyse
Folgende Daten sollen in einer einfachen Varianzanalyse untersucht
werden. Die Kennungen für die Gruppe müssen Zeichenketten sein. Wenn
nicht, muß der Vektor mit den Gruppenkennungen in einen factor
umgewandelt werden:
> gr <- factor(tf$gruppe)
Die Eingabe im Einzelnen:
wert gruppe
1 6 a
2 7 a
3 6 a
4 5 a
5 5 b
6 6 b
7 4 b
8 5 b
9 7 c
10 8 c
11 5 c
12 8 c
tafel <- read.table("1.dat")
tafel.aov <- aov(tafel$wert~tafel$gruppe)
> summary(tafel.aov)
Df Sum Sq Mean Sq F value Pr(>F)
tafel$gruppe 2 8.0000 4.0000 3.6 0.071 .
Residuals 9 10.0000 1.1111
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Zwei weitere Möglichkeiten, am Beispiel [7,Seite 388]:
Stichprobe 1 2 3
Werte 3 4 8
7 2 4
7 6
3
Eingabe der Daten
> wert <- c(3,7,4,2,7,3,8,4,6)
> gruppe <- c(1,1,2,2,2,2,3,3,3)
> gruppe<-factor(gruppe)
Variante 1:
> anova(lm(wert ~ gruppe))
Analysis of Variance Table
Response: wert
Df Sum Sq Mean Sq F value Pr(>F)
gruppe 2 6.8889 3.4444 0.6889 0.5379
Residuals 6 30.0000 5.0000
Variante 2:
> oneway.test(wert ~ gruppe,var.equal=T)
One-way analysis of means
data: wert and gruppe
F = 0.6889, num df = 2, denom df = 6, p-value = 0.5379
>
7.2 Bartlett Test
Der Bartlett-Test überprüft die Annahme, daß Streuungen in den Gruppen
annähernd gleich sind. Die Daten aus [9,Seite 222, 228] werden
in die Datei bartlett1.dat eingegeben:
Wert Gruppe
13 1
9 1
15 1
5 1
25 1
15 1
3 1
9 1
6 1
12 1
42 2
24 2
41 2
19 2
27 2
8 3
24 3
9 3
18 3
9 3
24 3
12 3
4 3
9 4
12 4
7 4
18 4
2 4
18 4
Dann erfolgt die Auswertung mit
> tafel <- read.table("bartlett1.dat",header=TRUE)
> werte <- tafel$Wert
> gruppen <- factor(tafel$Gruppe)
> bartlett.test(werte, gruppen)
Bartlett test for homogeneity of variances
data: werte and gruppen
Bartlett's K-squared = 1.6187, df = 3, p-value = 0.6552
>
Die Nullhypothese zur Gleichheit der Varianzen wird demnach nicht abgelehnt.
8 ,,Goodness of fit''-Tests (Anpassungstests an bekannte Verteilungen)
8.1 Kolmogorov-Smirnov Test für die Güte der Anpassung
Daten aus [10,S. 186]:
> da <- c(0.79, 0.68, 0.75, 0.73, 0.69, 0.77, 0.76, 0.74, 0.73, 0.68,
+ 0.72, 0.75, 0.71, 0.76, 0.69, 0.72, 0.70, 0.77, 0.71, 0.74)
> ks.test(da, "pnorm", mean=mean(da), sd=sqrt(var(da)))
One-sample Kolmogorov-Smirnov test
data: da
D = 0.0901, p-value = 0.997
alternative hypothesis: two.sided
Die Normalverteilungshypothese kann nicht verworfen werden.
8.2 Shapiro-Wilk-Test
Beispiel aus [11,Seite 398-399]: sind die Gewinne an
Gewicht von sieben Ratten (70, 85, 94, 101, 107, 118, 132)
normalverteilt?
> werte <- c(70, 85, 94, 101, 107, 118, 132)
> shapiro.test(werte)
Shapiro-Wilk normality test
data: werte
W = 0.998, p-value = 1
Es spricht nichts gegen die Annahme, daß die zugrundeliegende Verteilung
normalverteilt ist.
8.3 Überprüfung von Genotyp-Häufigkeiten auf Vorliegen eines
Hardy-Weinberg-Gleichgewichts
Beispiel aus [12]. Beobachtete Aufspaltungsziffern für Br(a/a),
Br(a/b), Br(b/b):
Br(a/a) | Br(a/b) | Br(b/b) |
1 | 22 | 86 |
> library(genetics)
...
> hpa5 <- c(rep("A/A",1),rep("A/B",22),rep("B/B",86))
> g1 <- genotype(hpa5)
> summary(g1)
Number persons typed: 109 (100%)
Allele Frequency: (2 alleles)
Count Proportion
A 24 0.11
B 194 0.89
Genotype Frequency:
Count Proportion
A/A 1 0.01
B/A 22 0.20
B/B 86 0.79
Heterozygosity (Hu) = 0.1977574
Poly. Inf. Content = 0.1767463
> HWE.chisq(g1,simulate.p.value=FALSE,correct=FALSE)
Pearson's Chi-squared test
data: tab
X-squared = 0.0985, df = 1, p-value = 0.7536
Warning message:
Chi-squared approximation may be incorrect in: chisq.test(tab, ...)
>
Das gleiche Ergebnis erhält man auch bei manualler Berechnung: Für ein
dialleles System bestimmt man die Häufigkeit der Allele (im Br(a/b)-Beispiel
ergibt sich 24/(2*109) für Br(a) und 194/(2*109) für Br(b)) und berechnet
die zu erwartenden Aufspaltungshäufigkeiten, wenn p die Häufigkeit des
Allels
a und q die Häufigkeit des Allels b ist, gemäß p2, 2pq, q2, die
entsprechenden Häufigkeiten sind 1,32110, 21,35780 und 86,32110. Dann
bestimmt man für alle drei Genotypen den Ausdruck [((B−E)2)/E] (B:
beobachtete Häufigkeit, E: erwartete Häufigkeit) und berechnet die Summe
dieser drei Werte, die für die Bewertung der Nullhypothese in diesem Beispiel
wie χ2 mit einem Freiheitsgrad verteilt ist
[7,Seite 251]. Im Beispiel ergibt das
χ2 = 0,09855. Damit gibt es keinen Anlaß für die Ablehnung der
Nullhypothese.
Die Verwendung
einer anderen Funktion aus dem genetics-Paket ergibt:
> HWE.exact(g1)
Exact Test for Hardy-Weinberg Equilibrium
data: g1
N11 = 86, N12 = 22, N22 = 1, N1 = 194, N2 = 24, p-value = 1
9 Vergleich von zwei unabhängigen Stichproben mit dem t-Test
Die Syntax lautet:
> t.test(x,y)
Oft befinden sich die Daten beider zu vergleichenden Gruppen in einem Vektor
und müssen aus diesem angand der Gruppenkennungen in einem zweiten Vektor
,,extrahiert'' werden. Beispiel aus P. Armitage und G. Berry: Statistical
Methods in Medical Research, S 113.
0 1.0
0 1.3
1 1.3
1 1.4
1 1.9
1 2.0
2 2.1
2 2.6
3 2.9
Dateneingabe: der Vektor werte enthält die Daten, der Vektor
gr definiert die Zugehärigkeit zu Gruppen 1, 2:
> werte <- c(0,0,1,1,1,1,2,2,3,1,1.3,1.3,1.4,1.9,2,2.1,2.6,2.9)
> gr <- c(rep(1,9), rep(2,9))
> gr <- factor(gr)
> t.test(werte[gr==1], werte[gr==2])
Welch Two Sample t-test
data: werte[gr == 1] and werte[gr == 2]
t = -1.5753, df = 13.844, p-value = 0.1378
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.4440276 0.2218053
sample estimates:
mean of x mean of y
1.222222 1.833333
>
10 Vergleich zweier abhängiger Stichproben im t-Test
Beide Wertereihen werden in Vektoren eingetragen:
> a <- c(4, 3.5, 4.1, 5.5, 4.6, 6, 5.1, 4.3)
> b <- c(3, 3, 3.8, 2.1, 4.9, 5.3, 3.1, 2.7)
> t.test(a,b, paired = TRUE)
Paired t-test
data: a and b
t = 2.798, df = 7, p-value = 0.0266
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.1781177 2.1218823
sample estimates:
mean of the differences
1.15
11 Wilcoxon-Test für unabhängige Stichproben
Beispiel Sachs 1984 (6. Aufl.) Seite 234:
> werte <- c(7,14,22,36,40,48,49,52,3,5,6,10,17,18,20,39)
> gruppe <- c(rep(1,8), rep(2,8))
> gruppe <- factor(gruppe)
> wilcox.test(werte[gruppe==1], werte[gruppe==2], alternative="greater")
Wilcoxon rank sum test
data: werte[gruppe == 1] and werte[gruppe == 2]
W = 53, p-value = 0.01406
alternative hypothesis: true mu is greater than 0
>
Alternativhypothese: Werte in Gruppe 1 sind größer.
12 Friedman-Test
groups in der Dokumentation sind ,,Behandlungen'', blocks
sind Individuen, Stichprobengruppen oder Blöcke. Beide
Eingabemöglichkeiten werden demonstriert. Beispiel aus Sachs:
Gruppe 1 2 3 4
Block
1 27 23 26 21
2 27 23 25 21
3 25 21 26 20
Die erste Variante ist günstig für die manuelle Eingabe
> tafel <- c(27, 27, 25, 23, 23, 21, 26, 25, 26, 21, 21, 20)
> mtafel <- matrix(tafel, nrow=3)
> friedman.test(mtafel)
Friedman rank sum test
data: mtafel
Friedman chi-squared = 8.2, df = 3, p-value = 0.04205
Diese Form der eingabe ist geeignet für das Auslesen aus einer Datei:
> werte <- c(27, 23, 26, 21, 27, 23, 25, 21, 25, 21, 26, 20)
> gruppe <- c(rep(c(1,2,3,4),3))
> gruppe <- factor(grupppe)
> block <- c(rep(1, 4), rep(2, 4), rep(3, 4))
> block <- factor(block)
> friedman.test(werte, gruppe, block)
Friedman rank sum test
data: werte, gruppe and block
Friedman chi-squared = 8.2, df = 3, p-value = 0.04205
13 H-Test (Kruskal-Wallis)
Beispiel aus Sachs (1984):
A B C D
12.1 18.3 12.7 7.3
14.8 49.6 25.1 1.9
15.3 10.1 47.0 5.8
11.4 35.6 16.3 10.1
10.8 26.2 30.4 9.4
8.9
> kwdaten <- c(12.1, 14.8, 15.3, 11.4, 10.8)
> kwdaten <- c(kwdaten, 18.3, 49.6, 10.1, 35.6, 26.2, 8.9)
> kwdaten <- c(kwdaten, 12.7, 25.1, 47.0, 16.3, 30.4)
> kwdaten <- c(kwdaten, 7.3, 1.9, 5.8, 10.1, 9.4)
> kwgruppe <- c(rep(1,5), rep(2,6), rep(3,5), rep(4,5))
> kruskal.test(kwdaten, kwgruppe)
Kruskal-Wallis rank sum test
data: kwdaten and kwgruppe
Kruskal-Wallis chi-squared = 11.5302, df = 3, p-value = 0.009179
alternative Form der Eingabe
> gra <- c(12.1, 14.8, 15.3, 11.4, 10.8)
> grb <- c(18.3, 49.6, 10.1, 35.6, 26.2, 8.9)
> grc <- c(12.7, 25.1, 47.0, 16.3, 30.4)
> grd <- c(7.3, 1.9, 5.8, 10.1, 9.4)
> res <- kruskal.test(list(gra, grb, grc, grd))
> res
Kruskal-Wallis rank sum test
data: list(gra, grb, grc, grd)
Kruskal-Wallis chi-squared = 11.5302, df = 3, p-value = 0.009179
>
14 Überlebenszeit-Analyse
14.1 Beispieldaten aus der R-Implementation benutzen
Der Aufruf der Daten des Datensatzes leukemia aus dem
survival-Paket ist in [13,Seite 238-9]
beschrieben.
> library(survival)
> data(leukemia)
> leukemia
time status x
1 9 1 Maintained
2 13 1 Maintained
3 13 0 Maintained
4 18 1 Maintained
5 23 1 Maintained
6 28 0 Maintained
7 31 1 Maintained
8 34 1 Maintained
9 45 0 Maintained
10 48 1 Maintained
11 161 0 Maintained
12 5 1 Nonmaintained
13 5 1 Nonmaintained
14 8 1 Nonmaintained
15 8 1 Nonmaintained
16 12 1 Nonmaintained
17 16 0 Nonmaintained
18 23 1 Nonmaintained
19 27 1 Nonmaintained
20 30 1 Nonmaintained
21 33 1 Nonmaintained
22 43 1 Nonmaintained
23 45 1 Nonmaintained
> leuk.fit <- survfit(Surv(time, status)~x, leukemia)
> summary(leuk.fit)
Call: survfit(formula = Surv(time, status) ~ x, data = leukemia)
x=Maintained
time n.risk n.event survival std.err lower 95% CI upper 95% CI
9 11 1 0.909 0.0867 0.7541 1.000
13 10 1 0.818 0.1163 0.6192 1.000
18 8 1 0.716 0.1397 0.4884 1.000
23 7 1 0.614 0.1526 0.3769 0.999
31 5 1 0.491 0.1642 0.2549 0.946
34 4 1 0.368 0.1627 0.1549 0.875
48 2 1 0.184 0.1535 0.0359 0.944
x=Nonmaintained
time n.risk n.event survival std.err lower 95% CI upper 95% CI
5 12 2 0.8333 0.1076 0.6470 1.000
8 10 2 0.6667 0.1361 0.4468 0.995
12 8 1 0.5833 0.1423 0.3616 0.941
23 6 1 0.4861 0.1481 0.2675 0.883
27 5 1 0.3889 0.1470 0.1854 0.816
30 4 1 0.2917 0.1387 0.1148 0.741
33 3 1 0.1944 0.1219 0.0569 0.664
43 2 1 0.0972 0.0919 0.0153 0.620
45 1 1 0.0000 NA NA NA
> plot(leuk.fit)
> survdiff(Surv(time, status)~x, data=leukemia, rho=0)
Call:
survdiff(formula = Surv(time, status) ~ x, data = leukemia, rho = 0)
N Observed Expected (O-E)^2/E (O-E)^2/V
x=Maintained 11 7 10.69 1.27 3.40
x=Nonmaintained 12 11 7.31 1.86 3.40
Chisq= 3.4 on 1 degrees of freedom, p= 0.0653
>
Mit library(Survival) wird das Survival-Paket geladen,
data(leukemia) lädt den Datensatz, leukemia zeigt die
Daten an. Mit Surv() wird ein ,,Überlebenszeit-Objekt''
erstellt, survfit() berechnet und plot() zeichnet eine
Überlebenszeitkurve (n. Kaplan Meier). Mit survdiff() wird die
Signifikanz der Unterschiede zwischen den Überlebenszeiten2 bestimmt. In der Datentabelle ist
unter time die Überlebenszeit beschrieben, mit status
wird der censoring status3 beschrieben, unter
x findet sich der Faktor, der die Gruppen definiert (,,maintenance
chemotherapy given?'').
14.2 Kaplan Meier-Kurve, Daten aus eigener Datentabelle einlesen
Nach Daten aus [14,S. 75] soll eine
Überlebenszeitkurve erstellt werden. Die geordnete Liste der Beobachtungen
(Zeit in Tagen) sei:
1*,3,4*,5,5,6*,7,7,7*,8*
wobei mit * die censored Überlebenszeiten (,,Ereignis
noch nicht eingetreten'') seien. Die Daten werden in die Datei
km.dat eingetragen:
time status
1 1 0
2 3 1
3 4 0
4 5 1
5 5 1
6 6 0
7 7 1
8 7 1
9 7 0
10 8 0
Die Auswertung und das Plotten der Kurve erfolgt dann mit:
> library(survival)
> tafel <- read.table("km.dat", header=T)
> attach(tafel)
> fit <- survfit(Surv(time, status) ~ 1, data=tafel)
> plot(fit)
Anmerkungen:
(1) Die rechte Seite der ,,Modellformel'' in survfit() für das Plotten einer Gruppe muss
,,~ 1'' lauten [15,Seite 67].
(2) Wenn
> attach(tafel)
nicht aufgerufen würde, müßte die Lebenzeitanalyse mit
> fit <- survfit(Surv(tafel$time, tafel$status), data=tafel)
veranlaßt werden.
15 Statistische Funktionen
15.1 Chiquadrat-Verteilung
> 1-pchisq(3.84, df=1)
[1] 0.05004352
> 1-pchisq(30.58, df=15)
[1] 0.009993624
15.2 F-Verteilung
Berechnung des p-Werts:
> 1-pf(9.55, df1=2, df2=3)
[1] 0.05001422
Berechnung des F-Werts (gegeben p, df1, df2):
> qf(0.025,df1=16,df2=26,lower.tail=F)
[1] 2.359684
15.3 t-Verteilung (Student)
> 1-pt(12.706,df=1) # fuer einseitigen Test
[1] 0.0250004
> 2*(1-pt(12.706,df=1)) # fuer zweiseitigen Test
[1] 0.0500008
> 2*(1-pt(2.787,df=25)) # fuer zweiseitigen Test
[1] 0.01001021
15.4 Standardnormalverteilung
> (1-pnorm(1.64485)) # fuer einseitigen Test
[1] 0.05000037
> 2*(1-pnorm(1.64485)) # fuer zweiseitigen Test
[1] 0.1000007
> 2*(1-pnorm(3.2905, mean=0, sd=1)) # fuer zweiseitigen Test
[1] 0.001000095
16 Bestimmung einer nowendigen Stichprobengröße
16.1 Vergleich zweier Mittelwerte
Beispiel 6.3 aus [11]: Vergleich der Lungenfunktion von zwei
Gruppen von Männern (mit dem forced expiratory volume), zuvor bekannt
eine Standardabweichung von 0.5, vorgegeben Signifikanzniveau von 0.05 für den
zweiseitigen Test und eine Power von 0.8, wie groß müssen die Stichproben sein
für die Erkennung einer Differenz von 0.25:
> power.t.test(delta=0.25, sd=0.5, sig.level=0.05, power=0.8,
+ alternative="two.sided" )
Two-sample t test power calculation
n = 63.76576
delta = 0.25
sd = 0.5
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group
Damit wird eine Sichprobengröße von 64 pro Gruppe benötigt.
16.2 Vergleich zweier relativer Häufigkeiten
Beispiel 6.5 aus [11]: Vergleich zweier Behandlungen: die
eine Behandlung führt zu einer Erfolgsrate von 25%. Wenn eine alternative
Behandlungsmethode zu einer Erfolgsrate vom 35% führt, wieviel Individuen
müssen die zu untersuchenden Gruppen bei einem Signifikanzniveau von 0.05 und
einer Power von 90% umfassen?
> power.prop.test(p1=0.25, p2=0.35, sig.level=0.05,
+ alternative="two.sided", power=0.9)
Two-sample comparison of proportions power calculation
n = 439.2309
p1 = 0.25
p2 = 0.35
sig.level = 0.05
power = 0.9
alternative = two.sided
NOTE: n is number in *each* group
also werden pro Gruppe 440 Individuen benötigt.
17 Metaanalyse
Im folgenden wird ein Beispiel aus [16] nachvollzogen. In einer
Metaanalyse werden Gruppen von Patienten mit perioperativer autologer
Blutspende daraufhin untersucht, wie wahrscheinlich die Gabe von Fremdblut
bei einer Operation sein wird. Hierzu wird das relative Risiko bestimmt.
Verglichen werden diese Wahrscheinlichkeiten mit Kontrollgruppen ohne
präoperative autologe Blutspende4.
Die Eingabe in die Datei blood-tr1.txt:
name n.trt n.ctrl tr.trt tr.ctrl
busch 239 236 66 133
elawad 45 15 3 14
hedstrom 38 40 7 29
heiss1993 58 62 20 37
heiss1997 29 27 11 13
hoynck 131 137 30 84
kajikawa 10 21 1 13
lorentz 16 15 2 10
Die Eingaben5:
> library(rmeta)
> d1 <- read.table("blood-tr1.txt",header=TRUE)
> d1
name n.trt n.ctrl tr.trt tr.ctrl
1 busch 239 236 66 133
2 elawad 45 15 3 14
3 hedstrom 38 40 7 29
4 heiss1993 58 62 20 37
5 heiss1997 29 27 11 13
6 hoynck 131 137 30 84
7 kajikawa 10 21 1 13
8 lorentz 16 15 2 10
> b <- meta.DSL(n.trt, n.ctrl, tr.trt, tr.ctrl, data=d1,statistic="RR",names=name)
> summary(b)
Random effects ( DerSimonian-Laird ) meta-analysis
Call: meta.DSL(ntrt = n.trt, nctrl = n.ctrl, ptrt = tr.trt, pctrl = tr.ctrl,
names = name, data = d1, statistic = "RR")
------------------------------------
RR (lower 95% upper)
busch 0.49 0.39 0.62
elawad 0.07 0.02 0.21
hedstrom 0.25 0.13 0.51
heiss1993 0.58 0.38 0.87
heiss1997 0.79 0.43 1.45
hoynck 0.37 0.27 0.53
kajikawa 0.16 0.02 1.07
lorentz 0.19 0.05 0.72
------------------------------------
SummaryRR= 0.38 95% CI ( 0.26,0.54 )
Test for heterogeneity: X^2( 7 ) = 22.39 ( p-value 0.0022 )
Estimated random effects variance: 0.14
> plot(b)
>
Resultat: in der Gruppe mit präoperativer autologer Spende ist die
Wahrscheinlichkeit 0.38 für eine Fremdbluttransfusion verglichen mit der
Kontrollgruppe, anders formuliert: die präoperative autologe Spende reduziert
das Risiko einer Fremdbluttransfusion um 62%.
18 Erstellung von Grafiken
Im folgenden Abschnit sind absichtlich die erzeugten Grafiken nicht
wiedergegeben, die entsprechenden Quellcodefragmente sollen interaktiv am
Propmpt eigegeben werden und man sieht dann, was passiert. Gleichzeitig sollte
man die Dokumentation öffnen z. B. für Barplot mit
help(barplot htmlhelp=T).
Allmählich wird man dann weitergehend probieren können.
18.1 Säulendiagramme
Es soll ein Säulendiagramm mit den folgenden Werten erstellt werden:
> werte <- c(0.03, 0.04, 0.69, 1.25, 0.08, 0.058, 1.1)
Es erscheinen senkrechte Säulen in Regenbogenfarben, Farbe wird kontrolliert
durch einen Vektor, der col zugewiesen wird
> farben <- rep(gray(0.7), 7)
> barplot(werte, col=farben)
nun sind alle Säulen hellgrau. Das Diagramm soll waagrechte Balken aufweisen:
> barplot(werte, col=farben, horiz=T)
Die Säulen sollen schmaler werden.
> barplot(werte, col=farben, horiz=T, space=3)
Darstellung gruppierter Säulen:
> werte <- c(0.03, 0.04, 0.69, 0.09, 1.25, 0.08, 0.58, 0.3)
> wertematrix <- matrix(werte, nrow=2)
> wertematrix
[,1] [,2] [,3] [,4]
[1,] 0.03 0.69 1.25 0.58
[2,] 0.04 0.09 0.08 0.30
> barplot(wertematrix, col=farben, horiz=F)
Das ergibt ,,gestapelte Säulen'', nebeneinander gruppiere Säulen:
> barplot(wertematrix, col=farben, horiz=F, beside=T)
als horizontales gruppiertes Diagramm:
> barplot(wertematrix, col=farben, horiz=T, beside=T)
18.2 Farben und andere Grafikparameter
Dokumentation zu Grafikparametern erhält man mit
> help(par)
Eine Liste von Farben für col wird mit
> colors()
gezeigt.
18.3 Boxplots
Annähernd normalverteite Zufallswerte für 2 Gruppen werden erzeugt mit:
> werte1 <- rnorm(20, mean=10, sd=2)
> werte2 <- rnorm(50, mean=16, sd=2.6)
> boxplot(werte1, werte2)
Farbe, Bereich der y-Achse adjustieren
> boxplot(werte1, werte2, col="wheat1", ylim=c(0,25))
Gruppenbezeichnungen einfügen:
> boxplot(werte1, werte2, col="wheat1", ylim=c(0,30),
+ names=c("group A", "group B"))
Körper der Boxen verschieden färben:
boxplot(werte1, werte2, col=c("wheat1","tomato2"),
+ ylim=c(0,30), names=c("group A", "group B"), horizontal=T)
Ein Boxplot mit logarithmisch unterteilter y-Achse wird mit
dem Zusatz log="y" gezeichnet:
e <- 2.7182818
werte1 <- rnorm(50, mean=6, sd=0.9)
werte1 <- e**werte1
werte2 <- rnorm(50, mean=7, sd=1.0)
werte2 <- e**werte2
boxplot(werte1, werte2, col="tomato1", log="y", names=c("GR 1", "GR 2"))
Oft wird es notwendig sein, die Boxen in einem Boxplot zu gruppieren.
Dazu kann man die Option at nehmen. Im folgenden ein weiteres
komplettes Beispiel:
werteA1 <- rnorm(20, mean=12, sd=3)
werteA2 <- rnorm(20, mean=14, sd=3.4)
werteB1 <- rnorm(18, mean=15, sd=3.2)
werteB2 <- rnorm(22, mean=18, sd=4)
werteC1 <- rnorm(24, mean=10, sd=2.2)
werteC2 <- rnorm(17, mean=9, sd=2)
boxplot(
werteA1,
werteA2,
werteB1,
werteB2,
werteC1,
werteC2,
ylim=c(0,25),
at = c(1.2,1.8,3.2,3.8,5.2,5.8),
boxwex=0.4,
ylab="microgram Ig/ml",
names=c("A1","A2", "B1", "B2",
"C1","C2"), las=2)
abline(v=2.5)
abline(v=4.5)
Um die Positionen der Säulen gleichmäßig zu verteilen, müßte man
at = c(1,2,3,4,5,6) eintragen. Mit boxwex skaliert man die
Breite der Boxen.
Mit zuätzlichen Linien (abline(...)) kann die Gruppierung deutlicher
gemacht werden.
18.4 Scatterplots
Im folgenden Beispiel wird ein Plot nach und nach um weitere Punkte, Legenden
ergänzt.
a <- c(1,2,2.5)
b <- c(2.6,4,6)
plot(a,b, xlim=c(0,3), ylim=c(0,8),
ylab="Y-Werte", xlab="X-Werte", pch=2, col="red")
c <- c(1.1,1.6,1.8,2)
d <- c(2.1,3,3.4,6)
points(c,d, pch=3, type="b", col="blue")
e <- c(1.3,1.7,1.9)
f <- c(2.1,2.7,3.4)
points(e,f, pch=4, type="b", col="magenta")
legend(2.5,1.2,plot=T,legend=c("Group a","Group c","Group e"),pch=c(2,3,4),
col=c("red","blue","magenta"))
18.5 Stripchart (eindimensionaler Scatterplot)
Während ein Boxplot ein abstrahiertes Bild der Verteilungsform liefert, ist es
manchmal auch erwünscht, Originalwerte verschiedener Gruppen direkt zu
plotten. Dazu kann man einen Stripchart verwenden.
A <- c(1, 3.4, 4, 4.2, 6, 12, 17)
B <- c(4, 8, 12, 13, 14.6, 18.9, 21.4, 14)
stripchart(list(A, B))
Gruppennamen werden vergeben und mit anderen Werten für ylim werden die beiden
Punktesäulen etwas in die Mitte geholt.
stripchart(list(A, B), group.names=c("group A", "group B"), ylim=c(0.5,2.5))
Um die Beschriftung der y-Achse horizontal zu setzen und wird die Anweisung
par(las=1) eingefügt und um den linken Rand zu verbreitern,
par(mar=c(5,6,5,2). Damit sieht der Quellcode für die Erstellung des
Diagramms so aus:
A <- c(1, 3.4, 4, 4.2, 6, 12, 17)
B <- c(4, 8, 12, 13, 14.6, 18.9, 21.4, 14)
stripchart(list(A, B))
par(las=1, mar=c(5,6,5,2))
stripchart(list(A, B), group.names=c("group A", "group B"), ylim=c(0.5,2.5))
Weitere Einzelheiten sind in der Dokumentation unter ?stripchart zu
erfragen.
Eine Alternative besteht in der Verwendung von stripplot aus dem
lattice-Paket:
library(lattice)
A <- c(1, 3.4, 4, 4.2, 6, 12, 17)
B <- c(4, 8, 12, 13, 14.6, 18.9, 21.4, 14)
trellis.device(color=FALSE)
gruppen <- factor(c(rep("A", length(A)), rep("B", length(B))))
stripplot(c(A,B)~gruppen, ylab="Wert", horizontal=FALSE)
18.6 Linien verbinden Werte von beliebig vielen Individuen,
die zwei oder mehr Zeitpunkten oder Wiederholungen entsprechen
18.6.1 Allgemeine Konstruktion mit plot()
par(adj = 0.5,
xaxt ="n")
plot(
# erste Linie
c(1,2), c(2.262,0.014),
type = "b",
xlim = c(0.6,2.4),
ylim = c(0,3),
ylab = "O. D.",
xlab = ""
)
text(1,2.9, vfont=c("sans serif", "plain"), cex=1.4, labels="before absorption")
text(2,2.9, vfont=c("sans serif", "plain"), cex=1.4, labels="after absorption")
# 2. und weitere Linien
lines(type="b", c(1,2), c(0.710,0.214))
lines(type="b", c(1,2), c(0.333,0.041))
lines(type="b", c(1,2), c(1.575,0.042))
lines(type="b", c(1,2), c(2.555,0.063))
lines(type="b", c(1,2), c(2.319,0.283))
lines(type="b", c(1,2), c(0.745,0.063))
lines(type="b", c(1,2), c(0.401,0.023))
lines(type="b", c(1,2), c(0.309,0.052))
lines(type="b", c(1,2), c(0.414,0.040))
lines(type="b", c(1,2), c(0.604,0.042))
lines(type="b", c(1,2), c(0.564,0.021))
lines(type="b", c(1,2), c(0.571,0.043))
lines(type="l", c(0.4,2.6),c(0.2,0.2))
Diese Art Darstellung setzt voraus, daß mit par(xaxt="n") die
Beschriftung der x-Axhse ausgeschltet wird. Dann wird die erste Linie im
Plot-Kommando gezeichnet, und dann folgen die weiteren Linien mit
lines()-Befehlen.
18.6.2 Konstruktion ähnlicher Diagramme mit interaction.plot()
Im einfachsten Fall ist die Syntax des Befehls:
interaction.plot(x.factor, trace.factor, response). Als Beispiel sei
angenommen, daß Werte (values, im Diagramm mit ,,Intensity''
beschriftet) von vier Individuen (subject) zu drei
verschiedenen Zeitpunkten (time: t0, t30, t90) untersucht wurden:
values = c(10,12,13,8,3,6,15,16,11,1,2,0.3)
time = rep(c("t0","t30","t90"),4)
subject= c(rep("no 1", 3),rep("no 2", 3),rep("no 3", 3),rep("no 4", 3))
interaction.plot(time,subject,values, ylab="Intensity")
Ein ähnliches Beispiel ist in [2,Seite 124] besprochen. Dieser
Diagrammtyp ist meist geeignet zur Darstellung von Daten, die mit der
Friedman-Rangvarianzanalyse untersucht werden.
19 Verschiedenes
19.1 Berechnung der Determinante einer Matrix
Im Beispiel werde für die Berechnung der Determinante:
die Funktion det aus dem Paket base.
> da <- c(0,-3,-8,1,4,3,0,-15,-3)
> matda <- matrix(da, nrow=3, byrow=T)
> det(matda)
[1] 111
19.2 Datumsberechnungen
Berechung des Abstands zwischen zwei Daten in Tagen. Beispiel: berechnet
werden soll der Abstand in Tagen zwischen 13.12.2000 und 20.4.2001:
> library(date)
> mdy.date(4,20,2001)-mdy.date(12,13,2000)
[1] 128
>
Die Differenz beträgt also 128 Tage.
Bestimmung des Datums n Tage vor oder nach einem gegebenen Daten: gesucht sei
das Datum eines Tages 100 Tage nach dem 16.10.2007:
> mdy.date(10,16,2007)+100
[1] 24Jan2008
>
Offenbar werden Daten intern als Ganzzahlen mit dem 1.1.1960 als Wert ,,0''
dargestellt:
> as.integer(mdy.date(1,1,1960))
[1] 0
> as.integer(mdy.date(1,2,1960))
[1] 1
>
20 Anhang
20.1 Hinweise für R unter Linux
20.1.1 Installation
Das Basisprogramm wird für einige Linux-Distributionen in fertig kompilierter
Form angeboten . Die Installation wird am Shell-Prompt mit
den Kommandos6
sudo apt-get install r-base
sudo apt-get install r-recommended
ausgelöst.
Zur Installation von Pakete besorgt man sich von der CRAN-Webseite die
entsprechenden Quellcodedateien (.tar.gz) und installiert am Shell-Prompt
(z.B. combinat_0.0-6.tar.gz)
R CMD INSTALL combinat_0.0-6.tar.gz
was funktioniert, wenn sich die notwendigen Entwicklungswerkzeuge
einschließlich gcc und Fortran-Compiler auf dem Linux-System befinden.
Außerdem sollten sich die Bibliotheken an der Stelle nach einer typischen
Installation befinden (Einzelheiten s. Installationsanleitung).
20.1.2 R konfigurieren
Das Verhalten von R kann durch den Befehl options()
beeinflußt werden, dessen Argument eine Liste von Zuweisungen im Format
,,variable=wert'' enthält.Will man beispielsweise einen anderen
Browser (SeaMonkey)
für die Anzeige der Hilfedateien definieren, so kann das (Beispiel für
Linux) unter Verwendung des nach chromium mit
options(browser="/usr/bin/chromium-browser")
,,von Hand'' geschehen, während R läuft. Wenn man dagegen eine
Konfigurationsdatei anlegen will, um dies vor dem Start des Programms
einzutragen, kann das in Rprofile.site im
etc-Unterverzeichnis geschehen. Weitere Optionen zu
Konfigurationsdateien sind im Abschnitt ,,Customizing the environment''
im Kapitel ,,Writing your own functions'' des
Begleitdokuments ,,An introduction to R'' beschrieben.
20.1.3 PDF-, PostScript- und EPS-Dateien von Abbildungen erstellen.
Unter Windows kann eine Grafikdatei mit dem Menü des Grafikfensters
gespeichert werden. Um eine EPS-Datei mit der Linux-Implementation
von R zu erstellen, wird die Grafik im Grafikfenster erzeugt mit einem dazu
geeigneten Befehl, dann am R-Prompt
dev.copy2eps()
eingegeben. Es wird dann eine Datei mit dem Namen Rplot.eps angelegt.
Um eine PostScript-Datei (Name: foo.ps) zu erzeugen, ist vor
dem Aufruf des Bildes
postscript("foo.ps")
einzugeben, das Bild mit dem geeigneten Aufruf zu erzeugen und anschließend mit
dev.off()
der PostScript-Ausgabeprozeß zu schließen. Ähnlich kann eine Grafik in eine
PDF-Datei (aus.pdf) geschrieben werden:
pdf("aus.pdf")
Befehle eingeben, um die PDF-Datei zu erzeugen, mit
dev.off()
den Ausgabeprozeß schließen.
20.2 Allgemeine Befehle
Der Befehl library() zeigt die Liste der verfügbaren Bibliotheken
an, mit library(MASS) wird die Library MASS geladen. Mit
data() wird die Liste der Datensätze angezeigt,
data(geyser) lädt und geyser zeigt die Daten zu
geyser.
Ermitteln des aktuellen Arbeitsverzeichnisses: getwd(), wechseln des
aktuellen Verzeichnisses: setwd("LW:/verzeichnis"). Mit
list.files() werden die Dateien im aktuellen Verzeichnis angezeigt.
Es empfiehlt sich bei der Arbeit mit R, den Quellcode im Texteditor
zu bearbeiten und zum Ausführen den abgespeicherten Quellcode auszuführen:
source("scatter1.r"). Wenn der auzuführende Code ausgedruckt werden
soll, geschieht das mit source("scatter1.r", echo=T). Wenn der
ausgedruckte Quellcode ohne Prompt erscheinen soll, geben Sie bitte
source("scatter1.r", echo=T,prompt.echo="") ein.
20.3 Einstellen von Optionen des Systems
Beim Start ist die Anzeige der ,,nackten'' Hilfedatein etwas unkomfortabel.
Die ,,Compiled-HTML-Hilfe'' kann als Vorgabe eingestellt werden mit:
options(htmlhelp=T)
Eine Option kann dann abgefragt werden mit z. B.:
getOption('chmhelp')
[1] TRUE
Wenn die Genauigkeit des Systems bei der Ausgabe von Resultaten geändert
werden soll, kann das mit
options(digits=20)
geschehen.
20.4 Textdarstellung von R-Objekten als Textdatei speichern
und wieder einlesen
Die Datei morley.tab werde eingelesen mit
mm <- read.table("morley.tab")
mit mm können dann die Daten angesehen werden. Der Vektor
a <- c(1, 2, 3)
und mm sollen in eine Textdatellung von R-Objekten
geschrieben werden:
dump(c("mm","a"), file = "mm.dat")
Mit
source("mm.dat")
werden dann beide Objekte wieder eingelesen.
20.5 Datendateien für R formatieren
Das folgende awk-Skript in Abschnitt 20.6.1 erlaubt die
Umwandlung einer Datendatei in das Format für R (es wurde für
gawk, die GNU-Implementation von awk geschrieben).
20.6 Quellcodes
20.6.1 ToR.awk
BEGIN {
i = 1
}
{
a = $0
if (i == 1)
{
printf(" %s\n", a)
}
else if (i > 1)
{
if (length(a) > 0)
{
printf("%-4i %s\n",i-1, a)
}
}
i++
}
Der Aufruf schreibt den Inhalt von input.txt formatiert in
out.dat.
gawk -f ToR.awk input.txt > out.dat
benutzt werden.
20.6.2 Perc.R
#
# perc(arr, pc.wert) returns the `pc.wert's percentile of the data
# in the vecctor `arr'
# Reference: Reed et al., Clinical Chemistry 1971:17(4);275-84
#
perc <- function(arr, pc.wert) {
neu <- arr
neu <- sort(neu)
n <- length(neu)
exakt.perz <- (n+1)*pc.wert/100
hilfswert <- floor(exakt.perz);
untergr <- neu[hilfswert];
obergr <- neu[hilfswert + 1];
if ((hilfswert < 1) || ((hilfswert+1) > n))
{
print("Percentile is out of range");
return(invisible())
}
perzentil <- untergr + (obergr - untergr)*(exakt.perz - floor(exakt.perz))
return(perzentil)
}
21 Versionen
04.03.2018: Abschnitt
1 eingefügt, Abschnitt
20.1 gekürzt und aktualisiert, Abschnitt
14:
R-Code zum zweiten Beispiel
aktualisiert.
References
- [1]
-
An Introduction to R.
http://cran.r-project.org/doc/manuals/R-intro.html.
- [2]
-
Dalgaard P.
Introductory Statistics with R.
Springer, New York, Berlin, Heidelberg, Hong Kong, London, Milan,
Paris, Tokyo, 2002.
- [3]
-
Maindonald J & Braun J, Hg.
Data analysis and graphics using R - an example-based approach.
Cambridge University press, Cambridge, 2003.
- [4]
-
Reed AH, Berry RJ, & Mason WB.
Influence of statistical method on the resulting estimate of normal
range.
Clinical Chemistry 17(4):274-284, 1971.
- [5]
-
Sachs L & Hedderich J.
Angewandte Statistik.
Springer, Dordrecht, Heidelberg, London, New York, 13. Aufl., 2009.
- [6]
-
Bortz J, Lienert GA, & Boehnke K, Hg.
Verteilungsfreie Methoden in der Biostatistik.
Springer-Verlag, Berlin, 1990.
- [7]
-
Sachs L.
Angewandte Statistik.
Springer, Berlin, Heidelberg, New York, Tokyo, 6. Aufl., 1984.
- [8]
-
Sachs L & Hedderich J, Hg.
Angewandte Statistik. Methodensammlung mit R.
Springer, Berlin, Heidelberg, New York, 12. Aufl., 2006.
- [9]
-
Storm R.
Wahrscheinlichkeitsrechnung, mathematische Statistik und statistische
Qualitätskontrolle.
Fachbuchverlag, Leipzig, Köln, 1995.
- [10]
-
Hartung J, Elpelt B, & Klösener KH.
Statistik.
R. Oldenbourg, München, Wien, 8. Aufl., 1991.
- [11]
-
Armitage P & Berry G.
Statistical methods in medical research.
Blackwell Scientific Publications, London, 3. Aufl., 1994.
- [12]
-
Kiefel V, Santoso S, & Mueller-Eckhardt C.
The Br(a)/Br(b) alloantigen system on human platelets.
Blood 73:2219-2223, 1989.
- [13]
-
Krause A & Olson M.
Statistics and computing.
Springer, New York, Berlin, Heidelberg, 2. Aufl., 2000.
- [14]
-
Matthews DE & Farewell VT.
Using and understanding medical statistics.
Karger, Basel, 2. Aufl., 1988.
- [15]
-
Matthews DE & Farewell VT.
Using and understanding medical statistics.
Karger, Basel, 5. Aufl., 2015.
- [16]
-
Henry DA, Carless PA, Moxey AJ, O'Connell D, Forgie MA, Wells PS, & Fergusson
D.
Pre-operative autologous donation for minimising perioperative
allogeneic blood transfusion.
Cochrane database of systematic reviews: CD003602, 2001.
Footnotes:
1frei verfügbar unter
http://cran.r-project.org/
2bei
rho=0 wird der Log-Rank oder Mantel-Haenszel Test, bei
rho=1 die PETO & PETO Modifikation des
GEHAN-WILCOXON-Tests verwendet
3üblicherweise: 0=noch
lebend, 1: verstorben (Ereignis eingetroffen)
4n.trt: ,,n treated'',
Anzahl der mit autologer Spende ,,behandelten'' Patienten, n.ctrl,
Anzahl der Patienten in der Kontrollgruppe, tr.trt: Anzahl der
Patienten mit Fremdbluttransfusionen in der Gruppe mit autologer
präoperativer Blutspende, tr.ctrl: Anzahl der Patienten mit
Fremdbluttransfusionen in der Kontrollgruppe
5meta.DSL für Random effects (DerSimonian-Laird)
meta-analysis
6hier beschrieben am Beispiel Ubuntu
File translated from
TEX
by
TTH,
version 4.12.
On 4 Mar 2018, 14:40.