# # Sitzung des Stochastik-Praktikums, WS 2013/2014, JGU Mainz # 27.1.2014 ## Thema: Lineares Modell mit normalverteiltem Rauschen ## ############################ ## Beispiel: Lineare Regression ## Die Siegzeiten im 800m-Lauf der Herren bei den Olympischen Spielen 1920-2008: laeufe <- read.table("800m.txt", header=TRUE) laeufe laeufe <- cbind(laeufe, data.frame(Zeit=60*laeufe$Min+laeufe$Sek)) laeufe attach(laeufe) plot(Jahr, Zeit, ylab="Zeit [s]", main="Siegzeiten im olympischen 800m-Lauf") modell <- lm(Zeit ~ Jahr) modell abline(modell) # Regressionsgerade einzeichnen modell$coeff a <- as.double(modell$coeff[1]); b <- as.double(modell$coeff[2]) a + b*2012 # das wäre die Vorhersage für 2012 gewesen # tatsächliche Siegzeit 2012: 1:40.91 ( David Rudisha) summary(modell) # mehr Informationen zur Anpassung abfragen modell$residuals # die Residuen plot(Jahr, modell$residuals) # anschauen hist(modell$residuals, prob=TRUE) detach(laeufe) ## ############################ ## Beispiele zur faktoriellen Varianzanalyse ## Beispiel: Blutgerinnungszeit von Ratten unter 4 verschiedenen ## Behandlungen ("A"-"D") wurde gemessen. ## Lade die Blutgerinnungsdaten (gerinnung.txt muss im Arbeitsverzeichnis sein): rat<-read.table("gerinnung.txt",header=TRUE) rat # Die Gruppenmittelwerte: mean(rat$bgz[rat$beh=="A"]); mean(rat$bgz[rat$beh=="B"]); mean(rat$bgz[rat$beh=="C"]); mean(rat$bgz[rat$beh=="D"]); # Globalmittel: mean(rat$bgz) # Anova-Tafel für die Blutgerinnungsdaten mit R: rat.aov <- aov(bgz~beh, data=rat) summary(rat.aov) # Zur Fisher-Verteilung: # qf(a, df1=k1,df2=k2) berechnet das a-Quantil der Fisher-Verteilung # mit k1 und k2 Freiheitsgraden, d.h. diejenige Zahl q mit der # Eigenschaft P(F <= q)=a, wo F Fisher(k1,k2)-verteilt. # Beispiel: # Wie muss man q wählen, damit P(F <= q)=0.95 für Fisher(6,63)-verteiltes F? qf(0.95,df1=6,df2=63) # pf(x, df1=k1,df2=k2) berechnet die Verteilungsfunktion der # Fisher(k1,k2)-Verteilung, d.h. # P(F <= x)=pf(x, df1=k1,df2=k2) für Fisher(k1,k2)-verteiltes F. # Beispiel: Wie wahrscheinlich ist es, dass eine # Fisher(3,20)-verteilte Zufallsgröße einen Wert <=13.57 annimmt? pf(13.57, df1=3, df2=20) # Demnach ist folgendes die Wahrscheinlichkeit, dass eine # Fisher(3,20)-verteilte Zufallsgröße einen Wert >13.57 annimmt: 1-pf(13.57, df1=3, df2=20) # Übrigens: R berechnet dasselbe auch via pf(13.57, df1=3, df2=20, lower.tail=FALSE) ## Beispiel: ## 7 Labore haben jeweils in 10 Proben den Gehalt eines Wirkstoffs gemessen, ## Daten aus ## R.D. Kirchhoefer, ## Semiautomated method for the analysis of chlorpheniramine ## maleate tablets: collaborative study, ## J. Assoc. Off. Anal. Chem. 62(6):1197-1201 (1979), ## zitiert nach John A. Rice, ## Mathematical statistics and data analysis, 2nd ed., Wadsworth, 1995 # # Lese Daten ein (Datei chlorpheniraminmaleat.txt muss im # Arbeitsverzeichnis sein): chlor <- read.table("chlorpheniraminmaleat.txt") chlor # Beachte: Hier ist es wichtig, durch den Befehl as.factor # die Labornummer zum "Faktor" zu machen (d.h. zu einem Gruppenbezeichner # ohne numerische Bedeutung, andernfalls interpretiert R Labor als Zahl; # dies war für gerinnung.dat nicht erforderlich, da R Buchstaben # ("A", "B", etc) automatisch als Faktoren interpretiert): str(chlor) # str(Objekt) (="structure") zeigt die interne Struktur eines R-Objekts chlor$Labor <- as.factor(chlor$Labor) str(chlor) # Visualisiere die Daten titel1<-c("7 verschiedene Labors haben jeweils 10 Messungen", "des Chlorpheniraminmaleat-Gehalts von", "Medikamentenproben vorgenommen:") stripchart(Gehalt~Labor, data=chlor, xlab="Labor", ylab="Gehalt an Chlorpheniraminmaleat [mg]", vertical=TRUE, method="stack",offset=0.6, main=titel1) boxplot(chlor$Gehalt~chlor$Labor, xlab="Labor", ylab="Gehalt an Chlorpheniraminmaleat [mg]",add=TRUE) # Bestimme empirische Mittelwerte und deren Standardfehler # für die Labors mwe <- numeric(7) sfe <- numeric(7) for (i in 1:7) { mwe[i] <- mean(chlor$Gehalt[chlor$Labor==i]) sfe[i] <- sd(chlor$Gehalt[chlor$Labor==i])/sqrt(10) } # Trage auch Mittelwerte +/- Standardfehler ein: titel2<-c("7 verschiedene Labors haben jeweils 10 Messungen", "des Chlorpheniraminmaleat-Gehalts von", "Medikamentenproben vorgenommen:", "Mittelwerte +/- Standardfehler") stripchart(Gehalt~Labor, data=chlor,xlab="Labor", ylab="Gehalt an Chlorpheniraminmaleat [mg]", vertical=TRUE, method="stack",offset=0.6, main=titel2) for (i in 1:7) { points(c(i-0.3,i+0.3),c(mwe[i],mwe[i]),type="l",col="red") points(c(i-0.3,i+0.3),c(mwe[i]-sfe[i],mwe[i]-sfe[i]),type="l",col="red",lty=2) points(c(i-0.3,i+0.3),c(mwe[i]+sfe[i],mwe[i]+sfe[i]),type="l",col="red",lty=2) } # Anova-Tafel für die Chlorpheniraminmaleat-Messungen mit R: chlor.aov <- aov(Gehalt~Labor,data=chlor) summary(chlor.aov) # Betrachten wir paarweise t-Tests der Labor-Messwerte M <- matrix(ncol=7,nrow=7,dimnames=list(c("Lab1","Lab2","Lab3","Lab4","Lab5","Lab6","Lab7"),c("Lab1","Lab2","Lab3","Lab4","Lab5","Lab6","Lab7"))) for (i in 1:6) { for (j in (i+1):7) { M[i,j] <- round(t.test(chlor$Gehalt[chlor$Labor==as.character(i)],chlor$Gehalt[chlor$Labor==as.character(j)])$p.value,5) } } M ## und dasselbe mit R: pairwise.t.test(chlor$Gehalt, chlor$Labor, method="none") # ohne Korrektur für multiples Testen pairwise.t.test(chlor$Gehalt, chlor$Labor, method="bonferroni") # mit Bonferroni-Korrektur pairwise.t.test(chlor$Gehalt, chlor$Labor, method="holm") # mit Bonferroni-Holm-Korrektur ## Tukeys "honest significant differences": TukeyHSD(chlor.aov) plot(TukeyHSD(chlor.aov)) ## ############################ ## Beispiel zur polynomiellen Regression: ## mittlere Temperatur im August in Karlsruhe nach Jahren a<-read.table("karlsruhe_august.txt", header=TRUE) a attach(a) Zeit<-Jahr/1000 b<-data.frame(Temp=Temp, Zeit=Zeit, Zeit2=Zeit^2, Zeit3=Zeit^3) m1<-lm(Temp~Zeit, data=b) summary(m1) m2<-lm(Temp~Zeit+Zeit2, data=b) summary(m2) m3<-lm(Temp~Zeit+Zeit2+Zeit3, data=b) summary(m3) plot(Zeit, Temp, main="Mittlere Augusttemperatur in Karlsruhe", xlab="Jahr/1000",ylab="Temp") points(Zeit,fitted.values(m1),type='l',lty=1) points(Zeit,fitted.values(m2),type='l',lty=2) points(Zeit,fitted.values(m3),type='l',lty=3) legend(1.87,24,legend=c("Grad 1","Grad 2","Grad 3"),lty=c(1,2,3)) ## Bsp.: Prüfen wir (von Hand), ob der lineare Term (1.064*Jahr) ## in m1 signifikant beitraegt: length(Zeit) V.stern <- sum(residuals(m1)^2)/(206-2) Var.ti <- mean(Zeit^2)-mean(Zeit)^2 Var.ti; V.stern 1.064/sqrt(V.stern/(206*Var.ti)) qt(0.975,df=204) # nicht signifikant auf dem 5%-Niveau pt(abs(1.064/sqrt(V.stern/(206*Var.ti))),df=204, lower.tail=FALSE) # dies waere der p-Wert detach(a)