R-Skript aus dem Seminar im Sommersemester 2007

Achtung!!! Dieses Skript gibt es auch als R-Syntax-Datei (genauer: Skript-Datei) zum Herunterladen: Syntax-Datei mit Kommentaren herunterladen

Grundlagen
    Vektoren, Daten, Objekte
    Einfache Berechnungen und Funktionen
    Wahrscheinlichkeitsfunktionen
    Daten einlesen
    Fehlende Werte
Deskriptive Statistiken und Grafiken
    Deskriptive Statistiken
    Summary für einen einfachen Überblick
    Einfache Grafiken
    Mehrere Grafiken gleichzeitig darstellen
Einfache Verfahren und Poweranalysen
Varianzanalysen
Regression
    Lineare Regression
    Binär-logistische Regression
Multivariate Verfahren
    Clusteranalyse
    Faktorenanalyse
Grafische Oberflächen

Vektoren, Daten, Objekte

Vektoren und Objekttypen

# R's Grundelemente sind Vektoren
# Im Sinne von Zusammenfassungen gleichartiger Komponenten
# von denen es im Wesentlichen vier Arten gibt:
# Zahlen, logische Werte, Texte und Faktoren
# Eine Methode, Vektoren zu erstellen ist der "c"-Befehl
zahl <- c(1,2,3,4,5)
zahl
logisch <- c(T,F,F,T,T)
logisch
text <- c("eins","zwei","drei","vier","R ist lustig")
text
faktor <- factor(c("männlich","weiblich","weiblich","weiblich","männlich"))
faktor
class(zahl)
class(logisch)
class(text)
class(faktor)

Wiederholungen und Sequenzen

# Oft braucht man Zahlenreihen oder bestimmte Wiederholungen von Zahlen
1:5
2:8
seq(8,25,2)
rep(c(1,2,4),c(10,5,3))

Mehrdimensionale Objekte

# Vektoren können wiederum zu mehrdimensionalen Objekten -
# Matrizen, Listen und "Data Frames" kombiniert werden.
# Dabei können nur Listen und Data Frames verschiedene
# Objektarten enthalten!!!
# In der folgenden Matrix sind alle Vektoren als
# Zeichenfolgen umgewandelt
ma1 <- matrix(c(zahl,logisch,text,faktor),5,4)
ma1
class(ma1)
liste <- list(zahl,logisch,text,faktor)
liste
class(liste)

Indizierung - Teilmengen

# Elemente von Vektoren, Matrizen, Listen, ... können mit
# der Element-Nummer in eckigen Klammern einzeln abgerufen werden.
# Die Element-Nummern werden bei der vollständigen Anzeige
# des Objekts angezeigt:
ma1
ma1[2,3]
ma1[1,]
ma1[,1]
liste
liste[[3]][5]
# außerdem können Elemente durch logische Objekte mit den selben
# Dimensionen angesteuert werden. Dies ist sinnvoll, um bestimmte
# Werte auszulesen und eventuelle weiterzuverarbeiten
# (siehe z.B. unten Behandlung fehlender Werte)
ma1
ma1 == "1" # logische Matrix mit denselben Dimensionen wie ma1
ma1[ma1=="1"]
# Data Frames sind das übliche Format für Datensätze,
# wie wir sie von SPSS gewohnt sind.
# Dabei entspricht jede Zeile einer Beobachtungseinheit / Versuchsperson
daten <- data.frame(zahl,logisch,text,faktor)
daten
attributes(daten)
# Data Frames können Variablennamen, -labels,
# nach Wunsch auch Zeilennamen enthalten
names(daten) <- c("VP-Nummer","Raucher","Irgendein Text","Geschlecht")
names(daten)
daten$Geschlecht
# Da Data Frames oft einen ganzen Datensatz enthalten,
# ist es nützlich, die enthaltenen Variablen fest in den
# Arbeitsbereich einzubinden
attach(daten)
Geschlecht
Raucher
# Wichtig: Änderungen an den Daten werden dann nicht
# ins Data Frame geschrieben!
# Um sie zu ändern muss die Variable im Data Frame
# mit "$" angesteuert werden
Raucher <- rep(F,5)
Raucher
daten$Raucher
ls()
rm(Raucher)
# Einbindung wieder aufheben
detach(daten)
# In R können alle gewohnten Rechenarten verwendet werden:
2+3
5/7
2^4
sqrt(25)
exp(2)
# auf Vektoren oder mehrdimensionale Objekte angewendet,
# wird die Berechnung mit jedem Element ausgeführt:
zahl
zahl + 2
# es können natürlich auch zwei Vektoren oder Objekte verwendet werden:
# bei verschiedener Länge wird das kürzere Objekt wiederholt
zahl * zahl
zahl1 <- c(zahl, zahl*3)
zahl1
zahl1 - 2 * zahl

Grundfunktionen

# Berechnung der Fakultät wird als Produkt umgesetzt
prod(1:5)
# Zufallsziehungen mit der Funktion sample
sample(1:50,5)
sort(sample(1:49,7))
sample(c("Erfolg","Misserfolg"),10,replace=T,prob=c(0.9,0.1))
# n über k
choose(10,3)
data.frame(0:10,choose(10,0:10))

Eigenschaften von Verteilungen

# R hat verschiedenste Verteilungen implementiert
help.search("distribution")
# Zu jeder Verteilung bietet R vier Funktionen, mit denen
# (1) die Dichte / Punktwarscheinlichkeit
# (2) die kumulierte Wahrscheinlichkeit / Verteilungsfunktion
# (3) die Quantile der Funktion ausgegeben werden und
# (4) Zufallszahlen aus der Verteilung erzeugt werden kÃ¶nnen.
# Beispiel Normalverteilung
?dnorm
curve(dnorm(x),from=-4,to=4)
curve(pnorm(x),from=-4,to=4)
qnorm(.5)
qnorm(.95)
rnorm(10,mean=180,sd=10)
print(rnorm(10,mean=180,sd=10), digits=3)
# Beispiel Binomialverteilung
plot(0:50,dbinom(0:50,size=50,prob=.33),type="h")

Arbeitsverzeichnis

# Das Arbeiten mit Dateien kann erheblich vereinfacht werden
# Wenn man das Arbeitsverzeichnis setzt, in dem sich Dateien
# befinden bzw. wohin sie geschrieben werden sollen
# !!! die "\" von Windows müssen durch "/" ersetzt werden !!!
getwd()
#eigenes Arbeitsverzeichnissetwd("H:/Hiwis/Christoph/R_Katalog")

Daten in Textform

# Dateien mit Komma-separierten Werten
d1 <- read.csv("testdaten.txt",colClasses=c("numeric","logical","character","factor"))
d1

Daten aus SPSS

# Paket zum Datenimport verschiedener Formate
# Unter anderem SPSS: foreign
library(foreign)
d2 <- read.spss("Termin2.sav",to.data.frame=T)
d2
# Da SPSS das Dateiformat immer wieder leicht ändert
# und damit Probleme entstehen können,
# besser mit "portable"-Dateien arbeiten:
d3 <- read.spss("termin2.por",to.data.frame=T)
d3
# R verwendet das Kürzel "NA" für fehlende Werte aller Objektarten
zahl
zahl[2] <- NA
zahl
logisch
logisch[3:4] <- NA
logisch

Behandlung fehlender Werte

# Standardmäßig ist R sehr konservativ mit fehlenden Werten
# und gibt bei Operationen mit fehlenden Werten wiederum
# "NA" zurück:
mean(zahl)
# verschieden Funktionen verändern die Behandlung fehlender Werte
# z.B. "na.exclude". Die vorgenommene Behandlung sowie die
# betroffenen Werte werden als Attribute gespeichert.
zahl.naex <- na.exclude(zahl)
zahl.naex
mean(zahl.naex)
# die fehlenden Werte können natürlich auch direkt in Funktionen
# ausgeschossen werden:
mean(na.exclude(zahl))

Fehlende Werte aus SPSS-Daten

# Oft wählt man in SPSS bestimmte Zahlen zur Kodierung fehlender Werte
# Diese übernimmt R beim Einlesen standardmäßig NICHT!!!
# man beachte die 9 als NA-Kodierung in dieser Datei:
d3
# Diese Werte müssen also noch von Hand auf NA gesetzt werden.
# Hier sieht man auch, dass mittels Index ausgewählte Teile
# eines Data-Frames wiederum mittels Index angesteuert werden können,
# prinzipiell unbegrenzt oft
d3
d3[4:length(d3)]
d3[4:length(d3)][d3[4:length(d3)]==9]
d3[4:length(d3)][d3[4:length(d3)]==9] <- NA
d3
# Mittelwert
mean(d3$ALTER)
# Standardabweicheung
sd(d3$ALTER)
# Varianz
var(d3$ALTER)
# Median
median(d3$ALTER)
# Quartile
quantile(d3$ALTER)
# Andere Quantile, hier z.B. Dezile
quantile(d3$ALTER,seq(0,1,.1))
# Achtung bei NAs:
mean(d3$ITEM1)
## Lösung
mean(d3$ITEM1, na.rm=T)
mean(na.exclude(d3$ITEM1))

Summary für einen einfachen Überblick

summary(d3$ALTER)
summary(d3)
# Man beachte den zu entdeckenden Fehler in den Altersangaben
# (Maximum 454???) und die Darstellung für Faktoren (Geschlecht)
# wenn sie wie hier korrekt definiert sind

Einfache Grafiken

# Boxplot
boxplot(d3$ALTER)
## Den Fehler im Alter auf NA setzen
d3$ALTER[d3$ALTER > 100]
d3$ALTER[d3$ALTER > 100] <- NA
d3$ALTER[d3$ALTER > 100]
boxplot(d3$ALTER)
# Q-Q-Plots
# Ist die Variable normalverteilt?
qqnorm(d3$ALTER)
# Balkendiagramm
barplot(summary(d3$SEX))
levels(d3$SEX)
levels(d3$SEX)[1] <- "männlich"
levels(d3$SEX)
barplot(summary(d3$SEX), col=c("royalblue3","wheat1"))
## Mit addierten Säulen, quer
barplot(matrix(c(1,2,3,2,4,6),3,2), horiz=T, col=rainbow(3))
## Mit gruppierten Säulen
barplot(matrix(c(1,2,3,2,4,6),3,2), beside=T, col=rainbow(6))
# Farben für Diagramme bestimmen:
colors()
?colors
# schnelle LÃ¶sung für beliebig viele verschiedene Farben:
?rainbow
pie(rep(1,5), col=rainbow(5,end=.7))
pie(rep(1,20), col=rainbow(20,alpha=.6,gamma=.6))
# Kuchendiagramm
pie(summary(d3$SEX))
# Nach Gruppen dargestellte Verteilungen - Stripchart
stripchart(d3$ALTER~d3$SEX)
# Histogramme
hist(d3$ALTER)
hist(d3$ALTER, breaks=3)
hist(d3$ALTER, breaks=c(0,10,20,25,60,70,80))
# Allgemeine Grafikfunktion: plot
# Beispiel: Kumulative Häufigkeit
plot(sort(d3$ALTER),1:length(sort(d3$ALTER))/length(sort(d3$ALTER)),type="s")

Mehrere Grafiken gleichzeitig darstellen

# Boxplots nach Gruppen
boxplot(d3$ALTER ~ d3$SEX)
# Allgemeine Methode für kombinierte Grafiken
# Grafikparameter "multiple frames columnwise"
par(mfcol=c(2,1))
hist(d3$ALTER[d3$SEX=="männlich"],col="royalblue3")
hist(d3$ALTER[d3$SEX=="weiblich"],col="wheat1")
par(mfcol=c(1,1))

Einfache Verfahren und Poweranalysen

# T-Test mit einer Stichprobe
mean(d3$ALTER)
t.test(d3$ALTER,mu=30)
# Selbe Hypothese, nicht-parametrischer Test
wilcox.test(d3$ALTER,mu=30)
# T-Test mit zwei Stichprobeh
names(d3)
## Werte in zwei verschiedenen Variablen
t.test(d3$ITEM1,d3$ITEM2)
## Berechnung als Messwiederholung
t.test(d3$ITEM1,d3$ITEM2, paired=T)
## Werte sind Gruppen derselben Variablen
t.test(d3$ALTER~d3$SEX)
# F-Test zum Vergleich der Varianzen:
var.test(d3$ITEM1,d3$ITEM2)
var.test(d3$ALTER~d3$SEX)
# Nicht-Parametrischer Test mit zwei Stichproben
wilcox.test(d3$ITEM1,d3$ITEM2)
wilcox.test(d3$ITEM1,d3$ITEM2, paired=T)
wilcox.test(d3$ALTER~d3$SEX)
# Poweranalyse
# R berechnet den nicht gegebenen Wert von
# n, delta, sd, Signifikanzniveau und Power
power.t.test(delta=.5,sd=2,sig.level=.01,power=.9)
power.t.test(,n=450,delta=.5,sd=2,sig.level=.01)
power.t.test(,n=450,delta=.5,sd=2,sig.level=.01, type="paired")
power.t.test(,n=450,delta=.5,sd=2,sig.level=.01, type="one.sample")

Varianzanalysen

library(ISwR)
data(red.cell.folate)
attach(red.cell.folate)
summary(red.cell.folate)
#einfaktorielle ANOVA mit "folate" als aV und "ventilation" als uV
anova(lm(folate~ventilation))
#Voraussetzung: uV muss kategorial vorliegen:
data(juul)
attach(juul)
anova(lm(igf1~tanner)) # wrong!!!
detach(juul)
#zuerst die uV als faktor definieren:
juul$tanner <- factor(juul$tanner,labels=c("I","II","III","IV","V"))
data(juul)
attach(juul)
summary(tanner)
#paarweise Vergleiche anhand der Regressionskoeffizienten
summary(lm(folate~ventilation))
#paarweise t-Tests (multiple testing)
pairwise.t.test(folate,ventilation,p.adj="bonferroni")
pairwise.t.test(folate,ventilation)
#bei Verletzung der Varianzhomogenität: oneway test nach Welch
bartlett.test(folate~ventilation) #Test auf Varianzhomogenität
oneway.test(folate~ventilation)
#paarweise Vergleiche ohne Voraussetzung der Varianzhomogenität
pairwise.t.test(folate,ventilation,pool.sd=F)
#Grafik
xbar <- tapply(folate,ventilation,mean)
s <- tapply(folate,ventilation,sd)
n <- tapply(folate,ventilation,length)
sem <- s/sqrt(n)
stripchart(folate~ventilation,"jitter",jit=0.05,pch=16,vert=T)
arrows(1:3,xbar+sem,1:3,xbar-sem,angle=90,code=3,length=0.1)
lines(1:3,xbar,pch=4,type="b",cex=2)
#Kruskal-Wallis test (nonparametrisches Testverfahren)
kruskal.test(folate~ventilation)
#zweifaktorielle ANOVA
data(heart.rate)
attach(heart.rate)
heart.rate
anova(lm(hr~subj+time))
#Grafik
interaction.plot(time,subj,hr)
interaction.plot(ordered(time),subj,hr) #berücksichtigt ungleiche Zeitabstände
#Friedman test (nonparametric Testverfahren)
friedman.test(hr~time|subj,data=heart.rate)

Regression

Lineare Regression

allgemeines

library(ISwR)
data(thuesen)
attach(thuesen)
edit(thuesen)

lineare Regression

lm(short.velocity~blood.glucose)
summary(lm(short.velocity~blood.glucose))

Diagramm

plot(blood.glucose,short.velocity)
abline(lm(short.velocity~blood.glucose))
abline(1.1,0.022)

Residuen und vorhergesagte Werte

lm.velo <- lm(short.velocity~blood.glucose)
fitted(lm.velo)
resid(lm.velo)

Anzeige mit fehlenden Werten

options(na.action=na.exclude)
lm.velo <- lm(short.velocity~blood.glucose)
fitted(lm.velo)

Diagramm mit Residuen

plot(blood.glucose,short.velocity)
abline(lm(short.velocity~blood.glucose))
segments(blood.glucose,fitted(lm.velo),blood.glucose,short.velocity)

Diagramm: Residuen aufgetragen gegen vorhergesagte Werte

plot(fitted(lm.velo),resid(lm.velo))

Normalverteilung der Residuen

qqnorm(resid(lm.velo))

Konfidenz-und Vorhersageintervalle

predict(lm.velo)
predict(lm.velo,int="c")
predict(lm.velo,int="p")

Diagramm mit Konfidenzintervallen

pred.frame <- data.frame(blood.glucose=4:20)
pp <- predict (lm.velo, int="p", newdata=pred.frame)
pc <- predict (lm.velo, int="c", newdata=pred.frame)
plot(blood.glucose,short.velocity, ylim=range(short.velocity,pp,na.rm=T))
pred.gluc <- pred.frame$blood.glucose
matlines(pred.gluc, pc, lty=c(1,2,2), col="black")
matlines(pred.gluc, pp, lty=c(1,3,3), col="black")

Pearson

cor(blood.glucose, short.velocity, use="complete.obs")
cor(thuesen,use="complete.obs")
cor.test(blood.glucose,short.velocity)

Nonparametrische Korrelationskoeffizienten

cor.test(blood.glucose,short.velocity,method="spearman")
cor.test(blood.glucose,short.velocity,method="kendall")
#Dateneinlesen
library(foreign)
arbeit <- read.spss("Arbeit.por", use.value.labels=FALSE, to.data.frame=T)
arbeit
attach(arbeit)
#Berechnung der logistischen Regression 1
glm(ARBFAE~ALTER+BLUTHQ+SUBJEKTI,binomial)
arbfae.logit <- glm(ARBFAE~ALTER+BLUTHQ+SUBJEKTI,binomial)
summary(arbfae.logit)
#Gruppenunterschiede in uVs
boxplot(ALTER~ARBFAE)
boxplot(SUBJEKTI~ARBFAE)
boxplot(BLUTHQ~ARBFAE)
#logistische Regression 2
arbfae2.logit <- glm(ARBFAE~ALTER+SUBJEKTI,binomial)
summary(arbfae2.logit)
#logistische Regression 3: nur Alter
arbfae3.logit <- glm(ARBFAE~ALTER,binomial)
#Plot der logistischen Kurve: Alter
a = arbfae3.logit$coefficients[1]
b = arbfae3.logit$coefficients[2]
curve(1/(1+exp(-1*(a+b*x))), from=25, to=60, col="blue", xlab="uV", ylab="Wahrscheinlichkeit", main="")
points(ALTER,ARBFAE, pch=16)
#Plot der logistischen Kurve: subjektives Wohlbefinden
arbfae4.logit <- glm(ARBFAE~SUBJEKTI,binomial)
a = arbfae4.logit$coefficients[1]
b = arbfae4.logit$coefficients[2]
curve(1/(1+exp(-1*(a+b*x))), from=10, to=40, col="blue", xlab="uV", ylab="Wahrscheinlichkeit", main="")
points(SUBJEKTI,ARBFAE, pch=16)
#Unterschied zwischen beobachteten und vorhergesagten Werten
fitted(glm(ARBFAE~ALTER))
#Einfluß der uVs
anova (arbfae.logit, test="Chisq")

Daten vorbereiten

d3
# Eine Untergruppe der Variablen auswählen
d4 <- d3[4:18]
d4
# Da nicht Zusammenhänge zwischen Versuchspersonen, sondern
# zwischen den Variablen untersucht werden, data.frame
# transponieren (Für die Clusteranalyse)
d4t <- t(d4)
d4t

Clusteranalyse

#Distanzmatrix berechnen
d4t.dist <- dist(d4t)
d4t.dist
#Clusteranalse mit der Distanzmatrix
d4t.clust <- hclust(d4t.dist)
d4t.clust
#Baumdiagramm der Clusteranalse
plot(d4t.clust)

Faktorenanalse

#NAs ausschließen
d4.na <- na.exclude(d4)
d4.na
nrow(d4)
nrow(d4.na)
#Hauptfaktoren:
d4.fact <- factanal(d4.na,factors=6)
d4.fact
summary(d4.fact)
#Hauptkomponenten
d4.pcomp <- princomp(d4.na)
d4.pcomp
summary(d4.pcomp)
plot(d4.pcomp)
biplot(d4.pcomp)
pairs(d4.pcomp$scores[,1:6], gap=0)
pairs(d4.pcomp$scores[,1:3])

Grafische Oberflächen

# Starten des R-Commanders
install.packages("Rcmdr", dependencies=TRUE)
library(Rcmdr)
# RKward gibt es momentan nur unter Linux.
# Auch sind nur wenige Befehle über die grafische Oberfläche installiert.

Sektionen