##*************************************************** ##* Ad-Oculos-Projekt-Video GLM, Generalized Linear ##* Models (verallgemeinerte lineare Modelle) ##* ##* Datum: 07.03.2021 ##* Version: 1.0 ##* ##* Projekt-Seite: https://www.faes.de/ad-oculos/ ##* Günter Faes, spv@faes.de ##*************************************************** ## Pakete laden: library(MASS) # ML-Schätzung, 2. Beispiel ## Beispieldaten laden: data(trees) ## Beispieldaten trees: # 31 Beobachtungen, 3 Merkmale: # Girth : numerisch Baumdurchmesser (eigentlich eher als Umfang) in Zoll # Height : numerisch Höhe in Fuß # Volume : numerisch Volumen des Holzes in Kubikfuß ############################################################# # 1. Beispiel: # Maximale Likelihood für Parameter einer Verteilung schätzen, # das Ganze zu Fuß! # # Beispieldatensatz: trees # Merkmal: Girth ############################################################# # Beispieldaten als Tabelle und Struktur darstellen: View(trees) str(trees) # Übersicht Verteilungsparameter für trees: summary(trees) sd(trees$Girth) # Histogramm zur Darstellung der Verteilung: hist(trees$Girth, freq = FALSE, main = "Verteilung des Merkmals Girth", ylab = "Dichte", xlab = "Merkmal Girth") # Hinzufügen der geschätzten Dichteverteilung: lines(density(trees$Girth), col = "darkblue", type ="l", lwd = 2) # ***** Schätzung der Log-Likelihood für 2 Modelle: ***** # Modell 1: Normalverteilung, Mittelwert = 15, Standardbweichung = 5 # Likelihood: L_Modell_1 <- dnorm(trees$Girth, mean = 15, sd = 5) # Log-Likelihood Modell 1: LL_Modell_1 <- log(L_Modell_1) # Achtung: log = ln # Summe der Loglikelihood als Kennzahl: LL_Modell_1_summe <- sum(LL_Modell_1) cat("LL Modell 1, MW = 15, sd = 5: ", LL_Modell_1_summe) # Modell 2: Normalverteilung, Mittelwert = 13, Standardbweichung = 3 # Likelihood: L_Modell_2 <- dnorm(trees$Girth, mean = 13, sd = 3) # Log-Likelihood Modell 2: LL_Modell_2 <- log(L_Modell_2) # Summe der Loglikelihood als Kennzahl: LL_Modell_2_summe <- sum(LL_Modell_2) cat("LL Modell 2, MW = 13, sd = 3: ", LL_Modell_2_summe) # Darstellen als DataFrame: Baum <- data.frame(Durchmesser = trees$Girth, L_Modell_1 = L_Modell_1, # Likelihood Modell 1 LL_Modell_1 = LL_Modell_1, # Log-Likelihood Modell 1 L_Modell_2 = L_Modell_2, # Likelihood Modell 2 LL_Modell_2 = LL_Modell_2 # Log-Likelihood Modell 2 ) View(Baum) ############################################################# # 2. Beispiel: # Durchführung einer MLE für das Merkmal dbh.cm unter der # Annhame, dass es einer Normalverteilung folgt. # # Beispieldatensatz: trees ############################################################# # Berechnung der ML-Schätzer (MLE) Mittelwert und Standardabweichung # fitdistr des MASS-Paketes: MLE_Girth <- fitdistr(trees$Girth, densfun = "normal") # Ausgabe Mittelwert (mean), Standradabweichung (sd) und der dazugehörigen # geschätzten Standardfehler (in Klammern): MLE_Girth # Ausgabe der MLE, über die die Verteilungsparameter geschätzt wurden: MLE_Girth$loglik ############################################################# # 1. GLM-Beispiel: # # Datensatz: trees # Y: Volume (Baumvolumen) # X: Girth (Bumdurchmesser) # Verteilungsannahme: Normalverteilt # # Modell: Y ~ X, family = gaussian ############################################################# # Überblick über den Zusammenhang zwischen Baumdurchmesser und Volumen: plot(Volume ~ Girth, data = trees, main = "Zusammenhang zwischen Baumdurchmesser und Volumen", xlab = "Baumdurchmesser", ylab = "Volumen", pch = 3, # Darzustellendes Symbol cex = 2, # Größe des Symbols col = "darkgreen") # Farbe des Symbols # Durchführung der GLM: M1 <- glm(Volume ~ Girth, family = gaussian, data = trees) # Ausgabe des Modells über summary: summary(M1) # Einzeichnen der Regressionsgeraden in obigen Plot: lines(trees$Girth, M1$fitted.values, lwd = 2, col = "blue") ## Wo finde ich fitted.values? Führen Sie folgende Anweisung aus und gehen durch die Ausgabe! # str(M1) ############################################################# # 2. GLM-Beispiel: # # Datensatz: Elefanten.csv # POOLE, J.H.: Mate guarding, reproductive success and female choice in African elephants. # In: Animal Behaviour 37 (1989), S. 842–849 # # Y: Paarung (Anzahl erfolgreicher Paarungen über einen Zeitraum von 8 Jahren) # X: Alter (... des Elefantens) # Verteilungsannahme: Poisson # # Modell: Y ~ X, family = poisson ############################################################# # Daten laden: Elefanten <- read.csv2("Elefanten.csv") View(Elefanten) # Überblick über den Zusammenhang zwischen Baumdurchmesser und Volumen: plot(Paarung ~ Alter, data = Elefanten, main = "Zusammenhang zwischen Elefantenalter und erfolgreichen Paarungen", xlab = "Alter", ylab = "Paarungserfolg", pch = 1, # Darzustellendes Symbol cex = 1, # Größe des Symbols col = "darkblue") # Farbe des Symbols # Durchführung der GLM: M2 <- glm(Paarung ~ Alter, family = poisson, data = Elefanten) # Ausgabe des Modells: summary(M2) # Einzeichnen der Anpassung in obigen Plot: lines(Elefanten$Alter, M2$fitted.values, lwd = 2, col = "red") # Wie hoch ist der Paarungserfolg, wenn der Elefant ein Alter von 43 Jahren erreicht hat? Alter <- c(43) a <- M2$coefficients[1]; a # Schnittpunkt b <- M2$coefficients[2]; b # Steigung Paarungserfolg <- exp(a + b * Alter); Paarungserfolg # --- Skript-Ende ---