################################################################################ ## Nichtparametrische Regression (GAMs) ## Über die nichtparametrische Regression wird eine flexibele Modellierung des ## Einflusses metrischer Kovariablen auf eine abhängige Variable angestrebt. ## Diese Methode kann zum Einsatz kommen, wenn eine Lineariserung kurvlinearer ## Zusammenhänge nicht ausreichend erscheint. ## ## Skript zum Video "Nichtparametrische Regression (GAMs)" ## ## Ad-Oculos-Projekt, https://www.faes.de/ad-oculos/ ## Günter Faes ## Version 1.0, 02.11.2022 ## R-Version: 4.2.1 ## RStudio-Version: RStudio 2022.07.1+554 "Spotted Wakerobin" ## ## Siehe auch die Playlist "Regression/Korrelation". ## ################################################################################ ############################################################################### # # Einstieg in das Thema # # Das folgende Beispiel zeigt die Anpassung über kubische Splines an die # beobachteten Daten. # Die Anpassung wird über die Funktionen spline und smooth.spline aus dem # Paket stats geschätzt. cars ist der Beispieldatensatz. # # Das Beispiel ist der R-Dokumentation zu smooth.spline entnommen. # ############################################################################### # Ausgabe der anzupassenden Daten: plot(dist ~ speed, data = cars, main = "Beispiele zur Spline-Anpassung", xlab = "Geschwindigkeit in mph", ylab = "Bremsweg in ft", sub = "Daten wurden in den 1920er Jahren erhoben!") # Eine rohe Anpassung über spline, Methode fmm (kubische Anpassung): cars.spl <- spline(cars$speed, cars$dist, method = "fmm") lines(cars.spl, col = "DarkMagenta", lty = 4) # Durchführung geglätteten Anpassung über smooth.spline: cars.sm_spl <- with(cars, smooth.spline(speed, dist)) cars.sm_spl # Grafische Darstellung der Anpassung: lines(cars.sm_spl, col = "blue") # Durchführung 2 erneuter Anpassungen mit geänderten Parametern: # Interessant ist hier der Parameter df (equivalent degrees of freedom) # als Vorgabe der Parameterzahl Basisfunktionen (Knoten). # df = 5 ss5 <- smooth.spline(cars[,"speed"], cars[,"dist"], df = 5) # df = 10 ss10 <- smooth.spline(cars[,"speed"], cars[,"dist"], df = 10) ss5; ss10 # Grafische Ausgabe der Anpassungen: lines(ss5, lty = 2, col = "red") lines(ss10, lty = 3, col = "darkgreen") legend(5,120,c(paste("Funktionsoptimiertes df =",round(cars.sm_spl$df,1)), "Anpassung mit df = 5", "Anpassung mit df = 10", "Demo: Rohe Anpassung über spline"), col = c("blue","red", "darkgreen", "DarkMagenta"), lty = 1:4, text.col = c("blue","red", "darkgreen", "DarkMagenta"), bg = "lightgray") ############################################################################### ## Pakete: library(mgcv) # Stellt Methoden der nichtparametrische Regression zur Verfügung library(gamair) # Beispieldaten zur GAM, Datensatz "brain" library(rgl) # 3D-Darstellung ############################################################################### # # 1. Beispiel # ############################################################################### ## Beispieldaten: # Für dieses Beispiel werden die Daten wie in der Doku zu gam des mgcv-Paketes # simuliert: set.seed(2) ## gamSim-Argumente: # eg : Spezifiziert, welche Beispieldaten simuliert werden sollen # (3: Beispiel mit kontinuierlichen Variablen) # n : Anzahl der simulierten Daten # dist : Legt den Verteilungstyp fest # scale: Rauschlevel SimDaten <- gamSim(eg = 3, n = 400, dist = "normal", scale = 1) View(SimDaten) # Grafischer Überblick: # Ausgabe der Funktion f(z) = 0.2 * x2^11 * (10 * (1 - x2))^6 + 10 * (10 * x2)^3 * (1 - x2)^10: plot(SimDaten$f, main = "Simulierte Beispieldaten und zugrundeliegender Funktionsverlauf", ylab = "Y-Beobachtung = f(z) * x1 + e", xlab = "Beobachtungs-Index", type = "p", col = "darkgrey") # Ausgabe der simulierten abhängigen Variable y <- f(z) * x1 + e: points(SimDaten$y, col = "black") ## Modellschätzung: Modell_1 <- gam(y ~ s(x2), method = "GCV.Cp", data = SimDaten) Modell_1b <- gam(y ~ s(x2, k=20), method = "GCV.Cp", data = SimDaten) summary(Modell_1) summary(Modell_1b) # AIC-Kriterium, ist ein oft genutztes Modellgütekriterium: Modell_1$aic Modell_1b$aic #Diagnostik, k-Index: gam.check(Modell_1) gam.check(Modell_1b) ## Darstellen des geschätzten Zusammenhanges: points(Modell_1$fitted.values, col = "blue") points(Modell_1b$fitted.values, col = "green") legend("topright", c("Grau: Abbildung f(z)", "Schwarz: Simulierte Beispieldaten y <- f(z) * x1 + e", "Blau: Geschätzter Zusammenhang Modell 1", "Grün: Geschätzter Zusammenhang Modell 1b"), text.col = c("darkgrey", "black", "blue", "green")) # Plot-Funktion des mgcv-Paketes zur Darstellung des Standardfehlers: plot.gam(Modell_1, main = "Geschätzter Funktionsverlauf", xlab = "Beobachtungs-Index", ylab = "Geschätzte Y-Werte", shade = TRUE) # Modellvorhersage mit neu generierten X2-Daten: set.seed(42) SimDaten_neu <- gamSim(eg = 3, n = 400, dist = "normal", scale = 1.5) View(SimDaten_neu) pred_neueD_Modell_1 <- predict.gam(Modell_1, SimDaten_neu[3], data = , type = "link", se = TRUE) se_u <- pred_neueD_Modell_1$fit - pred_neueD_Modell_1$se.fit se_o <- pred_neueD_Modell_1$fit + pred_neueD_Modell_1$se.fit # Dataframe zur Darstellung erstellen: pred_neueD_Modell_1_df <- data.frame(se_unten = se_u, Modellvorhersage = pred_neueD_Modell_1$fit, se_oben = se_o, SimDaten[3]) View(pred_neueD_Modell_1_df) # Grafische Darstellung: plot(pred_neueD_Modell_1_df$Modellvorhersage, main = "Modellvorhersage Modell 1, Simulation", sub = "Neue simulierte X-Werte wurden verwendet!", ylab = "Vorhergesagte Y-Werte", xlab = "Index der simulierten Werte", col = "blue") points(pred_neueD_Modell_1_df$se_unten, col = "darkgrey") # SE unten points(pred_neueD_Modell_1_df$se_oben, col = "darkgrey") # SE oben ## Ein wenig grafische Diagnostik: # Grafisch: par(mfrow = c(2,1)) # 2 Grafiken in einer Darstellung # Darstellung der Residuen als QQ-Plot: qq.gam(Modell_1, main = "1. Modell: Darstellung der Residuen als QQ-Plot") # Darstellung der Residuen: plot(Modell_1$residuals, main = "1. Modell: Darstellung der Residuen") par(mfrow = c(1,1)) # Normale grafische Darstellung ############################################################################### # # 2. Beispiel # ############################################################################### ## Beispieldaten "brain": # Ein Dataframe mit 5 Spalten und 1567 Zeilen. # Jede Zeile bezieht sich auf ein "Voxel" (Datenelement in einem 3-dimensionalen Gitter) # des Bildes. Die Spalten sind: # X: Voxelposition auf der horizontalen Achse. # Y: Voxelposition auf der vertikalen Achse. # medFPQ: Median von drei wiederholten "Fundamental Power Quotient"-Werten am Voxel, # dies ist das Hauptmaß für die Hirnaktivität. # Regionalcode: Gibt an, zu welcher von mehreren Hirnregionen das Voxel gehört. Die Regionen werden # von den Experimentatoren definiert. # 0: Basisregion # 1: Region von Interesse # 2: Region die durch den experimentellen Stimulus aktiviert wurde # NA: Voxel ohne Zuordnung # meanTheta: Mittlere Phasenverschiebung am Voxel, über drei Messungen. # # Auszug aus der Studien-Zusammenfassung: # Die funktionelle Magnetresonanztomographie (FMRI) misst die physiologische Reaktion des # menschlichen Gehirns auf experimentell kontrollierte Stimulation. In einem periodisch # angelegten Experiment ist es von Interesse, auf einen Unterschied im Timing (Phasenverschiebung) # der Reaktion zwischen zwei anatomisch unterschiedlichen Hirnregionen zu testen. data(brain) View(brain) brain <- brain[brain$medFPQ > 5e-3,] # 2 Extremwerte entfernen View(brain) # Grafischer Überblick der beobachteten Daten: plot3d(brain$X, brain$Y, brain$medFPQ, xlab = "X", ylab = "Y", zlab = "medFPQ", main = "Analyse des Beispieldatensatzes brain", col = "blue") ## Modellschätzung: # k: Basisdimension, siehe "choose.k". Modell_2a <- gam(medFPQ ~ s(Y, X, k = 100), method = "GCV.Cp", data = brain) # k = 100 Modell_2b <- gam(medFPQ ~ s(Y, X, k = 300), method = "GCV.Cp", data = brain) # k = 300 # # ----------- Ausgabe Modell 2a: ---------------- summary(Modell_2a) gam.check(Modell_2a) # # ----------- Ausgabe Modell 2b: ---------------- summary(Modell_2b) gam.check(Modell_2b) # Darstellen des geschätzten Zusammenhanges: plot3d(brain$X, brain$Y, Modell_2a$fitted.values, col = "green", add = TRUE) plot3d(brain$X, brain$Y, Modell_2b$fitted.values, col = "red", add = TRUE) # Die Dokumantation zu gamair beschreibt weitere Anpassungsmöglichkeiten # unter "ch7 Code for Chapter 7: GAMs in Practice: mgcv". ### Skript Ende