############################################################################### ## Inhalt: ## Fortsetzung von Kapitel 7.1 (Regressionsmodelle) des Kurses von ## Ruckdeschel & Kohl ############################################################################### ########################################################### ## NHANES-Studie ########################################################### ## Zur Illustration verwenden wir in dieser Sitzung die Daten der NHANES-Studie ## aus den Jahren 2005-2006. Die Daten können von der Homepage des Centers for ## Disease Control and Prevention (http://www.cdc.gov/nchs/nhanes.htm) ## heruntergeladen werden und liegen in SAS XPORT Format vor. Wir lesen diese ## Datei mit Hilfe der Funktion "read.xport" aus dem Paket "foreign" ein. library(foreign) ####################################### ## Wir beginnen mit den Untersuchungsdaten der Körpermaße Untersuchung1 <- read.xport(file = "nhanes2006Unt1.xpt") ## Wir kontrollieren, ob das Einlesen funktioniert hat und werfen zugleich einen ## ersten Blick auf die Daten. head(Untersuchung1) tail(Untersuchung1) str(Untersuchung1) ## Genauere Informationen zu den in diesem Data.frame enthaltenen Variablen ## finden sich unter: ## http://www.cdc.gov/nchs/data/nhanes/nhanes_05_06/bmx_d.pdf ## und ## http://www.cdc.gov/nchs/data/nhanes/nhanes_05_06/varexam_d.pdf names(Untersuchung1) ## Wir werden im Folgenden nur einen Teil der Variablen berücksichtigen. ## Diese sind: ## - SEQN = Identifikationsnummer (key-Variable) ## - BMXWT = Gewicht in kg ## - BMXHT = Größe in cm ## - BMXWAIST = Bauumfang in cm ## Wir reduzieren den Data.frame auf diese Variablen Untersuchung1 <- Untersuchung1[,c("SEQN", "BMXHT", "BMXWT", "BMXWAIST")] str(Untersuchung1) ####################################### ## Es folgen die Untersuchungsdaten des Blutdrucks Untersuchung2 <- read.xport(file = "nhanes2006Unt2.xpt") ## Wir kontrollieren, ob das Einlesen funktioniert hat und werfen zugleich einen ## ersten Blick auf die Daten. head(Untersuchung2) tail(Untersuchung2) str(Untersuchung2) ## Genauere Informationen zu den in diesem Data.frame enthaltenen Variablen ## finden sich unter: ## http://www.cdc.gov/nchs/data/nhanes/nhanes_05_06/bpx_d.pdf ## und ## http://www.cdc.gov/nchs/data/nhanes/nhanes_05_06/varexam_d.pdf names(Untersuchung2) ## Wir werden im Folgenden wieder nur einen Teil der Variablen berücksichtigen. ## Diese sind: ## - SEQN = Identifikationsnummer (key-Variable) ## - BPQ150A = gegessen in den letzten 30 Minuten ## - BPQ150B = Alkohol in den letzten 30 Minuten ## - BPQ150C = Kaffee in den letzten 30 Minuten ## - BPQ150D = Zigaretten in den letzten 30 Minuten ## - BPXPLS = 60s Puls ## - BPXSY1 = Systolischer Blutdruck (1. Messung) ## - BPXDI1 = Diastolischer Blutdruck (1. Messung) ## - BPXSY2 = Systolischer Blutdruck (2. Messung) ## - BPXDI2 = Diastolischer Blutdruck (2. Messung) ## - BPXSY3 = Systolischer Blutdruck (3. Messung) ## - BPXDI3 = Diastolischer Blutdruck (3. Messung) ## - BPXSY4 = Systolischer Blutdruck (optionale 4. Messung) ## - BPXDI4 = Diastolischer Blutdruck (optionale 4. Messung) ## Diese Variablen wurden bestimmt bei allen Personen im Alter 8 - 150 Jahre ## Wir reduzieren den Data.frame auf diese Variablen Untersuchung2 <- Untersuchung2[,c("SEQN", "BPQ150A", "BPQ150B", "BPQ150C", "BPQ150D", "BPXPLS", "BPXSY1", "BPXDI1", "BPXSY2", "BPXDI2", "BPXSY3", "BPXDI3", "BPXSY4", "BPXDI4")] str(Untersuchung2) ####################################### ## Es folgen die demographischen Daten Demographie <- read.xport(file = "nhanes2006Demo.xpt") ## Wir kontrollieren erneut das Einlesen. head(Demographie) tail(Demographie) str(Demographie) ## Genauere Informationen zu den in diesen Data.frames enthaltenen Variablen ## finden sich unter: ## http://www.cdc.gov/nchs/data/nhanes/nhanes_05_06/demo_d.pdf ## und ## http://www.cdc.gov/nchs/data/nhanes/nhanes_05_06/vardemo_d.pdf names(Demographie) ## Wir werden im Folgenden wieder nur einen Teil der Variablen berücksichtigen. ## Diese sind: ## - SEQN = Identifikationsnummer (key-Variable) ## - RIAGENDR = Geschlecht ## - RIDAGEYR = Alter zum Zeitpunkt der Untersuchung (in Jahren), Werte > 85 wurden auf 85 gesetzt! ## - RIDRETH1 = Rasse/Ethnische Gruppe ## - RIDEXPRG = liegt eine Schwangerschaft vor (Frauen im Alter von 8-59) ## - INDHHINC = geschätztes gesamtes Haushaltseinkommen Demographie <- Demographie[,c("SEQN", "RIAGENDR", "RIDAGEYR", "RIDRETH1", "RIDEXPRG", "INDHHINC")] str(Demographie) ####################################### ## Bei der Variable "SEQN" handelt es sich um eine Schlüsselvariable (key) ## Wir fassen mit Hilfe dieser Variable die 3 Data.frames zu einem Data.frame ## zusammen Daten <- merge(Untersuchung1, Untersuchung2, by = "SEQN", all = TRUE) Daten <- merge(Daten, Demographie, by = "SEQN", all = TRUE) ## Wir räumen auf rm(Untersuchung1, Untersuchung2, Demographie) dim(Daten) ## Es bleibt uns ein Datensatz mit 22 Variablen und 10348 Beobachtungen. ####################################### ## Wir wollen nun die Blutdruckwerte genauer untersuchen. ## Es gab bis zu 4 Messungen. Diese wollen wir in einem ersten Schritt ## zu einem Wert zusammenfassen. ## Wir beginnen mit einer Darstellung der Werte library(RColorBrewer) mypalette <- brewer.pal(4, "Dark2") boxplot(Daten[,10:17], col = rep(mypalette, each = 2), las = 2) ## Wir bilden Mittelwerte und Mediane für den systolischen und den diastolischen ## Blutdruck Sys.mean <- rowMeans(Daten[,c(10, 12, 14, 16)], na.rm = TRUE) Sys.med <- apply(Daten[,c(10, 12, 14, 16)], 1, median, na.rm = TRUE) Dia.mean <- rowMeans(Daten[,c(11, 13, 15, 17)], na.rm = TRUE) Dia.med <- apply(Daten[,c(11, 13, 15, 17)], 1, median, na.rm = TRUE) ## Wir vergleichen Mittelwert und Median plot(Sys.mean, Sys.med, pch = 20, main = "Systolischer Blutdruck", xlab = "Mittelwert", ylab = "Median") abline(a = 0, b = 1, lwd = 2, col = "orange") ## Es zeigen sich keine deutlichen Unterschiede plot(Dia.mean, Dia.med, pch = 20, main = "Diastolischer Blutdruck", xlab = "Mittelwert", ylab = "Median") abline(a = 0, b = 1, lwd = 2, col = "orange") ## Es zeigen sich deutlichere Auffälligkeiten als im Fall des systolischen Bludrucks ## Wir fügen die Mittelwerte der systolischen Blutdrucke zum Datensatz hinzu Daten <- data.frame(Daten, "BPXSYS" = Sys.mean) ####################################### ## Wir wollen nun den systolischen Blutdruck genauer Untersuchen. ## Wir erzeugen hierfür einen neuen Data.frame, bei dem die Einzelmessungen ## sowie die Variablen "SEQN" und "RIDEXPRG" unberücksichtigt bleiben. Daten1 <- Daten[,-c(1, 10:17, 21)] str(Daten1) ## Wir entfernen alle Beobachtungen, bei denen wenigstens eine der Variablen ## einen NA-Eintrag hat Daten1.na <- is.na(Daten1) ind1.na <- apply(Daten1.na, 1, any) Daten1 <- Daten1[!ind1.na, ] dim(Daten1) ## Es verbleibt ein Datensatz mit 13 Variablen und 7070 Beobachtungen ## Wir wandeln kategorielle bzw. ordinale Merkmale in Faktoren um. Daten1[,c(4:7, 9, 11, 12)] <- lapply(Daten1[,c(4:7, 9, 11, 12)], factor) str(Daten1) ####################################### ## Welche der im Datensatz enthaltenen Variablen sind Risikofaktoren für einen hohen ## Blutdruck? ########################################################### ## Achtung: ## 1. Schritt: Genauer Blick auf die Daten durch deskriptive Statistik ## (Summaries, Plots)!!! ## Nicht hier!!! ########################################################### ## Wir verwenden ein lineares Modell. ## Die Angabe der Formel "BPXSYS ~ ." bedeutet, dass wir alle Variablen im ## Datensatz "Daten1" berücksichtigen wollen. fit1 <- lm(BPXSYS ~ ., data = Daten1) summary(fit1) ## Einige Variablen sind signifikant andere nicht ... ## Wie finden wir das "beste" Modell? -> Modellwahl ## Lesen Sie dazu Abschnitt 7.1.1(g) im Skript von Ruckdeschel & Kohl ## DIAGNOSTIK!!! (Residuenplots, Test, etc.) - nicht hier! ####################################### ## Es gibt also eine Vielzahl von Kriterien, die zur Variablenauswahl ## verwendet (statistische Tests, R^2, AIC, BIC, etc.) werden können. Jedoch ## kann man nicht sagen, welches am besten ist! ## Generell gibt es drei Möglichkeiten ## - forward selection: Starte mit dem Nullmodell (nur Intercept) und füge ## schrittweise Variablen zum Modell hinzu solange das ausgewählte ## Kriterium erfüllt ist. ## - backward elimination: Starte mit dem vollen Modell und entferne schrittweise ## Variablen bis das das ausgewählte Modell erfüllt ist. ## - all subsets: Betrachte alle möglichen Teilmodelle. ####################################### ####################################### ## Bemerkung: ## 1. Da es unklar ist, wie die einzelnen Variablen zusammenhängen (Simpson-Paradoxon!) ## kombiniert man machmal auch "forward selection" und "backward selection" also ## bereits ausgewählte bzw. weggelassene Variablen können wieder wegfallen bzw. ## hinzugenommen werden. ## 2. Die Menge aller Teilmodelle erreicht schnell eine nicht mehr handhabbare ## Anzahl. ####################################### ####################################### ## Wir beginnen mit der "forward selection" und verwenden die Funktion "stepAIC" ## aus dem Paket "MASS" (ähnlich zur Funktion step aus Paket "stats") library(MASS) fit2 <- lm(BPXSYS ~ 1, data = Daten1) ## Nullmodell summary(fit2) fit3 <- stepAIC(fit2, scope = ~ BMXHT+BMXWT+BMXWAIST+BPQ150A+BPQ150B+BPQ150C+BPQ150D+BPXPLS+RIAGENDR+RIDAGEYR+RIDRETH1+INDHHINC, direction = "forward") summary(fit3) ## alles außer der Variable "BMXHT" wurde hinzugefügt ... ## DIAGNOSTIK!!! (Residuenplots, Test, etc.) - nicht hier! ####################################### ## Es folgt die "backward elemination" fit4 <- stepAIC(fit1, scope = ~ 1, direction = "backward") summary(fit4) ## lediglich die Variable "BMXHT" wird entfernt ## DIAGNOSTIK!!! (Residuenplots, Test, etc.) - nicht hier! ####################################### ## Kombination von forward selection und backward elemination fit5 <- stepAIC(fit1, scope = ~ 1, direction = "both") summary(fit5) ## erneut wird lediglich die Variable "BMXHT" entfernt ## DIAGNOSTIK!!! (Residuenplots, Test, etc.) - nicht hier! ####################################### ## Zur Modellwahl gibt es außerdem das Paket "leaps" -> alle Teilmodelle library(leaps) M <- model.matrix(BPXSYS ~ ., data = Daten1) fit6 <- leaps(x = M, y = Daten1$BPXSYS, strictly.compatible = FALSE) plot(fit6$size, fit6$Cp, pch = 20) abline(a = 0, b = 1, lwd = 2) ## Nochmal der gleiche Plot nun aber nur für fit6$size > 5 plot(fit6$size[fit6$size > 5], fit6$Cp[fit6$size > 5], pch = 20) abline(a = 0, b = 1, lwd = 2) ## Heuristik: Cp = Anzahl Variablen + 1 (hier: Cp = size) und möglichst ## wenig Variablen ## Hier also evtl. 6 oder auch 7 Variablen abs(fit6$Cp[fit6$size %in% c(7, 8)] - fit6$size[fit6$size %in% c(7, 8)]) auswahl <- which.min(abs(fit6$Cp[fit6$size %in% c(7, 8)] - fit6$size[fit6$size %in% c(7, 8)])) fit6$which[fit6$size %in% c(7, 8),][auswahl,] ########################################################### ## Wir werden im Folgenden nur die Variablen "BMXWT" (Gewicht in [kg]) ## und "BMWHT" (Größe in [cm]) verwenden. Daten2 <- Daten[,c("BMXWT", "BMXHT")] ## Wie wir oben bei der Ausgabe von "str" gesehen haben, enthält der Datensatz ## viele "NA" (not available) Einträge. Wir entfernen nun also alle Datensätze, ## bei denen wenigstens eine der beiden Variablen einen "NA" Eintrag enthält; ## d.h., unvollständig sind. Daten2.na <- is.na(Daten2) ind2.na <- apply(Daten2.na, 1, any) sum(ind2.na) Daten2 <- Daten2[!ind2.na,] dim(Daten2) ## Es verbleibt also ein Datensatz ("data.frame") mit 2 Variablen (Spalten) und ## 8949 Beobachtungen (Zeilen). ## Um uns die Arbeit zu erleichtern, fügen wir Daten zum Suchpfad hinzu. attach(Daten2) ####################################### ## Wir stellen nun in einem ersten Schritt die Daten in einem Scatterdiagramm ## graphisch dar. plot(BMXHT, BMXWT, pch = 20) ## Aufgrund der großen Anzahl von Daten ist es schwierig, die genauere Struktur ## der Datenwolke zu erkennen. Durch Verwendung der Funktion "smoothScatter" ## aus dem Bioconductor-Paket "geneplotter" kann man die Verteilung der ## Datenpunkte noch etwas anschaulicher machen. library(geneplotter) smoothScatter(BMXHT, BMXWT) ## Mittels der Parameter "nrpoints" und "colramp" kann man steuern, ob und wie ## die Datenpunkte eingezeichnet werden und welche Farbverläufe die Punktwolke ## zeigen soll; z.B. library(RColorBrewer) smoothScatter(BMXHT, BMXWT, nrpoints = 10, colramp = colorRampPalette(brewer.pal(9, "YlGnBu")), pch = 20) ####################################### ## Der Plot deutet bereits auf eine gewisse Heteroskedastizität hin. ## Wir berechnen die Regressionsgerade fit7 <- lm(BMXWT ~ BMXHT) summary(fit7) ## Wir fügen die Regressiongerade zum Scatterplot hinzu abline(fit7, lwd = 2) ## Der Residuenplot zeigt deutliche Strukturen (-> Heteroskedastizität!) plot(fit7, which = 1) ## Auch liegt eine deutliche Abweichung von der Normalverteilung vor! plot(fit7, which = 2) ####################################### ## Möglicher Ausweg: Datentransformation! ####################################### ####################################### ## Häufig stellt die Logarithmus-Transformation einen Ausweg dar smoothScatter(BMXHT, log(BMXWT), nrpoints = 10, colramp = colorRampPalette(brewer.pal(9, "YlGnBu")), pch = 20) ## Deutet immer noch auf eine Heteroskedastitzität hin ... ## Wir berechnen die Regressionsgerade fit8 <- lm(log(BMXWT) ~ BMXHT) summary(fit8) ## Wir fügen die Regressiongerade zum Scatterplot hinzu abline(fit8, lwd = 2) ## Der Residuenplot zeigt deutliche Strukturen (-> Heteroskedastizität!) plot(fit8, which = 1) ## Auch liegt eine deutliche Abweichung von der Normalverteilung vor! plot(fit8, which = 2) ####################################### ## Nächster Versuch: Andere Transformationen - Box-Cox-Potenztransformationen ## Hierfür verwenden wir die Funktion "boxcox" aus dem Paket "MASS" boxcox(BMXWT ~ BMXHT) boxcox(BMXWT ~ BMXHT, lambda = seq(-0.1, 0, 0.01)) ## Wir versuchen also eine Transformation mit lambda = -0.06 ## Frage: Wie ist dies bei der Darstellung zu interpretieren??? ## Interpretation wäre z.B. einfacher im Fall, dass lambda = 0 (-> Logarithmus) oder ## lambda = +/- 0.5 (-> Quadratwurzel) BMXWT.bc <- (BMXWT^(-0.06) - 1)/(-0.06) smoothScatter(BMXHT, BMXWT.bc, nrpoints = 10, colramp = colorRampPalette(brewer.pal(9, "YlGnBu")), pch = 20) ## Die Heteroskedastizität scheint kleiner geworden zu sein ... ## Wir berechnen erneut die Regressionsgerade fit9 <- lm(BMXWT.bc ~ BMXHT) summary(fit9) abline(fit9, lwd = 2) ## Der Residuenplot zeigt jedoch immer noch deutliche Strukturen (-> Heteroskedastizität!) plot(fit9, which = 1) ## Auch liegt nach wie vor eine deutliche Abweichung von der Normalverteilung vor! plot(fit9, which = 2) ####################################### ## Möglichkeiten: ## 1. Bilde Subgruppen ## Zum Beispiel: BMXHT <= 140 und BMXHT > 140 bzw. noch mehr Gruppen ## und führe für jede Subgruppe eine eigene Analyse durch. ## 2. Die Variablen erklären den Zusammenhang unzureichend, es gibt ## noch unbekannte Einflussfaktoren. ####################################### ########################################################### ## Ein spezielles lineares Modell liegt vor, falls alle Regressoren Faktoren ## sind. Man spricht dann von ANOVA (Analysis of Variance); siehe Abschnitt 7.1.1(c) ## im Skript von Ruckdeschel & Kohl. ########################################################### ####################################### ## Wir betrachten die Blutdruckdaten ## 1-Weg ANOVA plot(BPXSYS ~ RIAGENDR) ## -> Boxplots ## 1. Möglichkeit fit10 <- lm(BPXSYS ~ RIAGENDR, data = Daten1) anova(fit10) ## 2. Möglichkeit fit11 <- aov(BPXSYS ~ RIAGENDR, data = Daten1) summary(fit10) ## 3. Möglichkeit oneway.test(BPXSYS ~ RIAGENDR, data = Daten1, var.equal = TRUE) ## aber auch möglich, verschiedene Varianz (nicht mit "lm" bzw. "aov"!) oneway.test(BPXSYS ~ RIAGENDR, data = Daten1, var.equal = FALSE) ## Zusammenfassung der Ergebnisse model.tables(fit11, type = "means") ####################################### ## 2-Weg ANOVA par(mfrow = c(1, 2)) plot(BPXSYS ~ RIAGENDR + RIDRETH1) ## -> Boxplots ## 1. Möglichkeit fit11 <- lm(BPXSYS ~ RIAGENDR + RIDRETH1, data = Daten1) anova(fit11) ## 2. Möglichkeit fit11 <- aov(BPXSYS ~ RIAGENDR + RIDRETH1, data = Daten1) summary(fit11) ## Zwei Faktoren können aber auch interagieren. Die entsprechenden Modelle ## in R sind. fit12 <- aov(BPXSYS ~ RIAGENDR + RIDRETH1 + RIAGENDR:RIDRETH1, data = Daten1) summary(fit12) ## keine signifikante Interaktion! ## identisch zu fit13 <- aov(BPXSYS ~ RIAGENDR*RIDRETH1, data = Daten1) summary(fit13) ## keine signifikante Interaktion! ## Zusammenfassung der Ergebnisse model.tables(fit11, type = "means") ## Graphische Darstellung (Interaktionsplot) ## Schneiden sich die Linien deutet dies auf eine Interaktion hin. interaction.plot(Daten1$RIAGENDR, Daten1$RIDRETH1, Daten1$BPXSYS) ## bzw. interaction.plot(Daten1$RIDRETH1, Daten1$RIAGENDR, Daten1$BPXSYS) ## hier keine Interaktion zu erkennen ## oder einfacher library(HH) interaction2wt(BPXSYS ~ RIAGENDR*RIDRETH1, data = Daten1) ########################################################### ## Allgemeiner führt dies zu einer k-Weg ANOVA. ## Meist liegt die Situation vor, dass man speziell geplante Versuche durchführt, ## um den Einfluss der Faktoren festzustellen. ## Bsp: Neue Produktionsmethode, neue Maschine, neues Medikament, etc. ## Dies führt ins Gebiet der Versuchsplanung/Design of Experiments. ###########################################################