# Stochastik-Praktikum, WS 2019/2020, JGU Mainz # 20.1.2017 # ## ############################################################################## ## ## Thema: ## [Etwas zur] Lineare[n] Regression ## ######################################################### # Lineare Regression, ein Beispiel: # Rs cats-Datensatz library("MASS") data(cats) ?cats attach(cats) koerper <- Bwt[Sex=="M"]; herz <- Hwt[Sex=="M"] # (Wir betrachten nur die 97 Kater) plot(koerper, herz, xlab='Körpergewicht [kg]',ylab='Herzgewicht [g]') # wir passen die Ausgleichsgerade # Herzgewicht = a + b * Koerpergewicht # an: modell <- lm(herz ~ koerper) modell # die Koeffizienten: modell$coefficients # die Residuen: modell$residuals # zeichnen wir die Regressionsgerade ein: abline(modell, col='red') # (aequivalent ist hier:) abline(modell$coefficients[1], modell$coefficients[2]) # hist(modell$residuals) qqnorm(modell$residuals) # Praediktion: a <- modell$coefficients[1]; b <- modell$coefficients[2] neues.kg <- 3.75 predicted.hg <- b * neues.kg + a predicted.hg points(neues.kg, predicted.hg, pch=22, col='blue') # Bem.: Extrapolation kann natuerlich auch in der Anwendungssituation # "sinnlose" Werte ergeben, z.B. neues.kg <- 0.25 # ein sehr leichter Kater predicted.hg <- b * neues.kg + a predicted.hg neuer.h.wert <- 25 # (z.B. ein sehr "gechillter" Geier) predicted.m.wert <- bb[1]+bb[2]*neuer.h.wert predicted.m.wert # (das waere eine negatives vorhergesagtes Herzgewicht, # was natuerlich keinen Sinn hat ...) detach(cats) ## ###################################################### # # Woher kommt der Name "Regression"? # Francis Galtons Beobachtung der "Rueckkehr zum Mittelwert" ## Beispiel: K. Pearsons Daten über Größen von Vätern und Söhnen pearson <- read.table("pearson.dat", header=T) pearson[1:10,] attach(pearson) length(fheight) plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)") # Mittelwerte einzeichnen mean(fheight) # 67.6871 mean(sheight) # 68.68407 abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) mo <- lm(sheight ~ fheight) mo # und die Regeressionsgerade: abline(mo, lwd=2, col="red") # (empirische) Korrelation cor(fheight, sheight) # 0.5013383 sq.f <- sum((fheight-mean(fheight))^2)/length(fheight) sq.f # 7.527313 sqrt(sq.f) # 2.743595 sq.s <- sum((sheight-mean(sheight))^2)/length(sheight) sq.s # 7.915196 sqrt(sq.s) # 2.813396 cov.fs <- sum((fheight-mean(fheight))*(sheight-mean(sheight)))/length(fheight) cov.fs # 3.869739 cov.fs/(sqrt(sq.f)*sqrt(sq.s)) # Korrelationskoeff: 0.5013383 # uebrigens: die Daten (zumindest die Marginalien) sind nicht unplausibel # durch die Normalverteilung beschrieben: hist(fheight) qqnorm(fheight) hist(sheight) qqnorm(sheight) ## ######################################### ## grosse und kleine Vaeter anschauen: ## grosse: gr1 <- 71.5; gr2 <- 72.5 ausw.gr <- (fheight >= gr1) & (fheight <= gr2) sum(ausw.gr) farbe.gr <- "blue" plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)", type="n") abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) abline(0, 1) # Diagonale einzeichnen abline(mo, lwd=2, col="red") points(fheight[!ausw.gr], sheight[!ausw.gr], pch=20, cex=0.5) points(fheight[ausw.gr], sheight[ausw.gr], pch=20, cex=0.5, col=farbe.gr) points((gr1+gr2)/2, mean(sheight[ausw.gr]), pch=19, cex=1.3, col=farbe.gr) ##points((gr1+gr2)/2, mean(sheight[ausw.gr]), pch=19, cex=0.7, col="white") ## points(fheight, sheight, pch=20, cex=0.5) # Histogramme vergleichen kgr <- 57:80 kgr2 <- seq(from=58, to=80, by=2) hi1 <- hist(sheight, breaks=kgr, plot=F) hi2 <- hist(sheight[ausw.gr], breaks=kgr2, plot=F) plot(hi1$mids, hi1$density, xlim=c(55, 80), ylim=c(0,1.1*max(hi1$density, hi2$density)), xlab="", ylab="Dichte", type="l") abline(h=0, lty=1, col="darkgray") points(hi2$mids, hi2$density, type="l", col=farbe.gr) abline(v=mean(sheight), lty=2) abline(v=mean(sheight[ausw.gr]), lty=2, col=farbe.gr) ## zusaetzlich mit Rand-Histogrammen: ## ## gr-Vaetern werden markiert: x11(width=9) layout(matrix(c(1,2),nrow=1), widths=c(3,1), heights=c(1,1)) # Scatter: par(mar=c(4.1,4.1,4.1,1)) plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)", type="n") abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) abline(0, 1) # Diagonale einzeichnen abline(mo, lwd=2, col="red") points(fheight[!ausw.gr], sheight[!ausw.gr], pch=20, cex=0.5) points(fheight[ausw.gr], sheight[ausw.gr], pch=20, cex=0.5, col=farbe.gr) points((gr1+gr2)/2, mean(sheight[ausw.gr]), pch=19, cex=1.3, col=farbe.gr) abline(v=gr1, lty=3, col=farbe.gr); abline(v=gr2, lty=3, col=farbe.gr) lines(c((gr1+gr2)/2,82),rep(mean(sheight[ausw.gr]),2),col=farbe.gr,lty=2) # dann Histogramm: par(mar=c(4.1,0,4.1,0.2)) plot(hi1$density, hi1$mids, ylim=c(55, 80), xlim=c(0,1.1*max(hi1$density, hi2$density)), ylab="", xlab="", #xlab="Dichte", main=c("Dichte","(Gr. d. Sohns)"), type="l") abline(v=0, lty=1, col="darkgray") points(hi2$density, hi2$mids, type="l", col=farbe.gr) abline(h=mean(sheight), lty=2) abline(h=mean(sheight[ausw.gr]), lty=2, col=farbe.gr) ## ############################################### ## kleine: kl1 <- 64.5; kl2 <- 65.5 ausw.kl <- (fheight >= kl1) & (fheight <= kl2) sum(ausw.kl) farbe.kl <- "darkgreen" plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)", type="n") abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) abline(0, 1) # Diagonale einzeichnen abline(mo, lwd=2, col="red") points(fheight[!ausw.kl], sheight[!ausw.kl], pch=20, cex=0.5) points(fheight[ausw.kl], sheight[ausw.kl], pch=20, cex=0.5, col=farbe.kl) points((kl1+kl2)/2, mean(sheight[ausw.kl]), pch=19, cex=1.3, col=farbe.kl) ##points((kl1+kl2)/2, mean(sheight[ausw.kl]), pch=19, cex=0.7, col="white") ## points(fheight, sheight, pch=20, cex=0.5) # Histogramme vergleichen kkl <- 57:80 kkl2 <- seq(from=58, to=80, by=2) hi1 <- hist(sheight, breaks=kkl, plot=F) hi3 <- hist(sheight[ausw.kl], breaks=kkl2, plot=F) plot(hi1$mids, hi1$density, xlim=c(55, 80), ylim=c(0,1.1*max(hi1$density, hi3$density)), xlab="", ylab="Dichte", type="l") abline(h=0, lty=1, col="darkgray") points(hi3$mids, hi3$density, type="l", col=farbe.kl) abline(v=mean(sheight), lty=2) abline(v=mean(sheight[ausw.kl]), lty=2, col=farbe.kl) ## mit Rand-Histogrammen, # kl-Vaetern markiert: x11(width=9) layout(matrix(c(1,2),nrow=1), widths=c(3,1), heights=c(1,1)) # Scatter: par(mar=c(4.1,4.1,4.1,1)) plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)", type="p") abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) abline(0, 1) # Diagonale einzeichnen abline(mo, lwd=2, col="red") points(fheight[!ausw.kl], sheight[!ausw.kl], pch=20, cex=0.5) points(fheight[ausw.kl], sheight[ausw.kl], pch=20, cex=0.5, col=farbe.kl) points((kl1+kl2)/2, mean(sheight[ausw.kl]), pch=19, cex=1.3, col=farbe.kl) abline(v=kl1, lty=3, col=farbe.kl); abline(v=kl2, lty=3, col=farbe.kl) lines(c((kl1+kl2)/2,82),rep(mean(sheight[ausw.kl]),2),col=farbe.kl,lty=2) # dann Histogramm: par(mar=c(4.1,0,4.1,0.2)) plot(hi1$density, hi1$mids, ylim=c(55, 80), xlim=c(0,1.1*max(hi1$density, hi3$density)), ylab="", xlab="", #xlab="Dichte", main=c("Dichte","(Gr. d. Sohns)"), type="l") abline(v=0, lty=1, col="darkgray") points(hi3$density, hi3$mids, type="l", col=farbe.kl) abline(h=mean(sheight), lty=2) abline(h=mean(sheight[ausw.kl]), lty=2, col=farbe.kl) ## ############################################### ## alles zusammen: kl und gr x11(width=9) layout(matrix(c(1,2),nrow=1), widths=c(3,1), heights=c(1,1)) # Scatter: par(mar=c(4.1,4.1,4.1,1)) plot(fheight, sheight, main="1078 Größen von Vater und Sohn", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, xlab="x[Vater] : Größe des Vaters (in Inches)", ylab="x[Sohn] : Größe des Sohns (in Inches)", type="n") abline(h=mean(sheight), lty=2); abline(v=mean(fheight), lty=2) abline(0, 1) # Diagonale einzeichnen abline(mo, lwd=2, col="red") points(fheight[!(ausw.gr | ausw.kl)], sheight[!(ausw.gr | ausw.kl)], pch=20, cex=0.5) points(fheight[ausw.gr], sheight[ausw.gr], pch=20, cex=0.5, col=farbe.gr) points((gr1+gr2)/2, mean(sheight[ausw.gr]), pch=19, cex=1.3, col=farbe.gr) abline(v=gr1, lty=3, col=farbe.gr); abline(v=gr2, lty=3, col=farbe.gr) lines(c((gr1+gr2)/2,82),rep(mean(sheight[ausw.gr]),2),col=farbe.gr,lty=2) points(fheight[ausw.kl], sheight[ausw.kl], pch=20, cex=0.5, col=farbe.kl) points((kl1+kl2)/2, mean(sheight[ausw.kl]), pch=19, cex=1.3, col=farbe.kl) abline(v=kl1, lty=3, col=farbe.kl); abline(v=kl2, lty=3, col=farbe.kl) lines(c((kl1+kl2)/2,82),rep(mean(sheight[ausw.kl]),2),col=farbe.kl,lty=2) # dann Histogramme: par(mar=c(4.1,0,4.1,0.2)) plot(hi1$density, hi1$mids, ylim=c(55, 80), xlim=c(0,1.1*max(hi1$density, hi2$density)), ylab="", xlab="", #xlab="Dichte", main=c("Dichte","(Gr. d. Sohns)"), type="l") abline(v=0, lty=1, col="darkgray") points(hi2$density, hi2$mids, type="l", col=farbe.gr) abline(h=mean(sheight), lty=2) abline(h=mean(sheight[ausw.gr]), lty=2, col=farbe.gr) points(hi3$density, hi3$mids, type="l", col=farbe.kl) abline(h=mean(sheight[ausw.kl]), lty=2, col=farbe.kl) ## ################################# ## Andersherum: Vater als Funktion des Sohns mo.andersherum <- lm(fheight ~ sheight) mo.andersherum plot(sheight, fheight, main="1078 Größen von Sohn und Vater", xlim=c(55, 80), ylim=c(55, 80), pch=20, cex=0.5, ylab="x[Vater] : Größe des Vaters (in Inches)", xlab="x[Sohn] : Größe des Sohns (in Inches)") abline(v=mean(sheight), lty=2); abline(h=mean(fheight), lty=2) abline(mo.andersherum, lwd=2, col="violet") abline(0, 1) # Diagonale einzeichnen # detach(pearson) ## ###################################################### # baeume <- read.table("trees.csv", header=TRUE, sep=',') attach(baeume) pairs(baeume) plot(DBH_cm, Height_m) plot(log(DBH_cm), log(Height_m)) plot(DBH_cm, Crown_radius_m) plot(log(DBH_cm), log(Crown_radius_m))