###############################################################################
## Inhalt:
## Abschnitt 5.2 (Schätzen) des R-Kurses von Ruckdeschel & Kohl
###############################################################################

## Wir laden den Datensatz
Daten <- read.csv2(file = 'PatientenDaten.csv')
attach(Daten)

## Patientendaten mit folgenden Variablen
## Patient: Patientennummer (1 - 31)
## Gruppe: Ausprägung einer bestimmten Erkrankung (schwach, mittel, stark)
## CRP: C-reaktives Protein (Entzündungsmarker) Werte in [0,Inf)
## SIRS-Kriterien: (Systemic Inflammatory Response Syndrom) 0, 1, 2, 3, 4

## Zusätzlich führen wir die Variable logCRP ein.
logCRP <- log(CRP)


###############################################################################
## Lesen Sie bitte den Abschnitt 5.2.1 im Skript!
###############################################################################

###########################################################
## Momentenschätzer
###########################################################
## mean, sd, var 
## sowie skewness, kurtosis (vgl. e1071)

## Wir betrachten die Variable logCRP
mean(logCRP, na.rm = TRUE) ## empirisches erstes Moment = empirischer Erwartungswert
mean(logCRP^2, na.rm = TRUE) ## empirisches zweites Moment
mean(logCRP^3, na.rm = TRUE) ## empirisches drittes Moment
## usw.

## Oder: Zentrale Momente (zentriert am empirischen Erwartungswert)
## Die folgende Funktion berechnet das k-te zentrale Moment
zentralMoment <- function(x, k){
  mean((x-mean(x))^k)
}
zentralMoment(logCRP[!is.na(logCRP)], 2) ## empirische Varianz


###########################################################
## Maximum-Likelihood (ML) Schätzer
###########################################################
## Im Fall des ML-Schätzers ist als eine Maximierungsaufgabe zu lösen.
## vgl. auch Abschnitt 5.2.2 im Skript!
## 
## 1. Möglichkeit: Analytisch über die Ableitung
## 2. Möglichkeit: Numerisch mittels 'optimize' bzw. 'optim'
## Ein einfaches Beispiel
fun <- function(x) -(x-3)^2
optimize(fun, interval = c(0, 5), maximum = TRUE)
## oder
fun1 <- function(x) -(x[1] - 3)^2 - (x[2] - 1)^2
optim(par = c(1, 1), fn = fun1, control = list(fnscale = -1)) 
## Per default sucht 'optim' nach einem Minimum.
## Durch Angabe von 'fnscale = -1' erreicht man, dass nach einem Maximum
## gesucht wird.

## 3. Möglichkeit: Numerische Lösung als Nullstelle der Ableitung mittels 'uniroot'
fun2 <- function(x) (x+3)^3
uniroot(fun2, interval = c(-5, 0))

## Die Berechnung von ML-Schätzer für eine ganze Reihe von Verteilungen
## ist in der Funktion 'fitdistr' aus dem Paket 'MASS' implementiert
library(MASS)

#######################################
## 1. Beispiel: Normalverteilung
#######################################
## Wir betrachten die Variable logCRP
(res <- fitdistr(logCRP[!is.na(logCRP)], densfun = 'normal'))

## Im Fall der Normalverteilung sind die Maximum-Likelihood-Schätzer
## gerade das arithmetische Mittel (mean) und die Standardabweichung (sd)

## Wir vergleichen die nichtparametrische Dichteschätzung mit der
## Schätzung, die wir mit Hilfe der ML-Methode erhalten
curve(dnorm(x, mean = res$estimate[1], sd = res$estimate[2]), from = -1, to = 7, 
      n = 101, col = 'orange', ylab = 'Dichte', lwd = 2,
      main = 'Parametrische vs. nichtparametrische Dichteschätzung')
lines(density(logCRP[!is.na(logCRP)]), lwd = 2)


#######################################
## 2. Beispiel: Gammaverteilung
#######################################
## Wir betrachten die Variable CRP
(res1 <- fitdistr(CRP[!is.na(CRP)], densfun = 'gamma'))

## Wir vergleichen die nichtparametrische Dichteschätzung mit der
## Schätzung, die wir mit Hilfe der ML-Methode erhalten
curve(dgamma(x, shape = res1$estimate[1], rate = res1$estimate[2]), from = 0, to = 250, 
      n = 501, col = 'orange', ylab = 'Dichte', lwd = 2,
      main = 'Parametrische vs. nichtparametrische Dichteschätzung')
lines(density(CRP[!is.na(CRP)], from = 0), lwd = 2)


#######################################
## Außerdem sind in 'fitdistr' implementiert:
## Beta-, Cauchy-, chi^2-, Exponential-, F-, Geometrisch, Lognormal-, Logistisch, 
## Negativ Binomial, Poisson, t- und Weibull-Verteilung
#######################################


###########################################################
## Lesen Sie Abschnitt 5.2.3 (b)-(e) aus dem Skript!!!
###########################################################


## Ein sehr allgemeines Schätz-Prinzip ist auch das Minimum-Distanz (MD) Prinzip
## Dabei wird ein geeigneter Abstand zwischen Verteilungen minimiert.
## Einer der bekanntesten Schätzer ist

###########################################################
## Kolmogorov(-Smirnov) Minimum Distanz (KS) Schätzer
## KS-Schätzer sind in gewisser Hinsicht robust gegenüber Abweichungen 
## von den Modellannahmen
###########################################################
## Ist implementiert in der Funktion 'ksEstimator' im Paket ROptEst 
library(ROptEst)

#######################################
## 1. Beispiel: Normalverteilung
#######################################
## Wir betrachten die Variable 'logCRP'
(res2 <- ksEstimator(logCRP, distribution = Norm()))

## Wir vergleichen die nichtparametrische Dichteschätzung mit der
## Schätzung, die wir mit Hilfe des KS-Schätzers erhalten
curve(dnorm(x, mean = res2$mean, sd = res2$sd), from = -1, to = 7, 
      n = 101, col = 'orange', ylab = 'Dichte', lwd = 2,
      main = 'Parametrische vs. nichtparametrische Dichteschätzung')
lines(density(logCRP[!is.na(logCRP)]), lwd = 2)


#######################################
## 2. Beispiel: Gammaverteilung
#######################################
## Wir betrachten die Variable 'CRP'
(res3 <- ksEstimator(CRP, distribution = Gammad()))

## Wir vergleichen die nichtparametrische Dichteschätzung mit der
## Schätzung, die wir mit Hilfe des KS-Schätzers erhalten
curve(dgamma(x, shape = res3$shape, scale = res3$scale), from = 0, to = 250, 
      n = 501, col = 'orange', ylab = 'Dichte', lwd = 2,
      main = 'Parametrische vs. nichtparametrische Dichteschätzung')
lines(density(CRP[!is.na(CRP)], from = 0), lwd = 2)


###########################################################
## M-Schätzer 
## M-Schätzer sind in gewisser Hinsicht (optimal-)robust gegenüber Abweichungen 
## von den Modellannahmen
###########################################################
## Im Fall der Normalverteilung gibt es Implementation des sog. Huber-Schätzers
## in 'MASS' ('huber' bzw. 'hubers') und 'robustbase' ('huberM')
library(MASS)
huber(logCRP)
(res4 <- hubers(logCRP))

library(robustbase)
huberM(logCRP)

#######################################
## 1. Beispiel: Normalverteilung
#######################################
## Wir betrachten die Variable 'logCRP'
## Wir vergleichen die nichtparametrische Dichteschätzung mit der
## Schätzung, die wir mit Hilfe der M-Methode erhalten
curve(dnorm(x, mean = res4$mu, sd = res4$s), from = -1, to = 7, 
      n = 101, col = 'orange', ylab = 'Dichte', lwd = 2,
      main = 'Parametrische vs. nichtparametrische Dichteschätzung')
lines(density(logCRP[!is.na(logCRP)]), lwd = 2)


#######################################
## 2. Beispiel: Gammaverteilung
#######################################
## Mir ist keine verwendbare Implementation in R bekannt. Es gibt jedoch 
## Arbeiten von Alfio Marazzi dazu ...

## Generell funktioniert der M-Schätzer Ansatz nur gut in einfachen Modellen. 
## Hingegen ist der folgende Ansatz in beliebigen (glatt-parametrisierten) 
## Modellen möglich!


###########################################################
## Asymptotisch Lineare (AL) Schätzer
## AL-Schätzer sind in gewisser hinsicht robust gegenüber Abweichungen von den 
## Modellannahmen
###########################################################
## Sind in den Paketen ROptEst, RobLox, ROptRegTS und RobRex implementiert.
## Es fehlt noch an einem einfacheren User-Interface, auch müssen die
## Algorithmen noch hinsichtlich der Rechenzeit optimiert werden.
library(RobLox)

#######################################
## 1. Beispiel: Normalverteilung
#######################################
## Wir betrachten die Variable logCRP
(res5 <- roblox(logCRP[!is.na(logCRP)], initial.est = 'med', returnIC = TRUE))

## Wir vergleichen die nichtparametrische Dichteschätzung mit der
## Schätzung, die wir mit Hilfe der AL-Schätzer erhalten
curve(dnorm(x, mean = res5$mean, sd = res5$sd), from = -1, to = 7, 
      n = 101, col = 'orange', ylab = 'Dichte', lwd = 2,
      main = 'Parametrische vs. nichtparametrische Dichteschätzung')
lines(density(logCRP[!is.na(logCRP)]), lwd = 2)

## Die sog. Influenzkurve, mißt den Einfluss einer einzelnen Beobachtung
## auf den Schätzer
plot(res3$optIC)

#######################################
## 2. Beispiel: Gammaverteilung
#######################################
## Generell benötigt man einen Startschätzer für die Berechnung der 
## AL-Schätzer
(est0 <- ksEstimator(x = CRP, Gammad()))

## Anlegen des 'robusten Modells'
(RobG <- InfRobModel(center=GammaFamily(scale = est0$scale, shape = est0$shape), 
                neighbor=ContNeighborhood(radius=0.25)))

## Berechnen der Influenzkurve
IC <- optIC(model=RobG, risk=asMSE())
plot(IC)

## Berechnung des Schätzers
(est1 <- oneStepEstimator(CRP[!is.na(CRP)], IC=IC, start=est0))

## Wir vergleichen die nichtparametrische Dichteschätzung mit der
## Schätzung, die wir mit Hilfe des KS-Schätzers erhalten
curve(dgamma(x, shape = est1[2], scale = est1[1]), from = 0, to = 250, 
      n = 501, col = 'orange', ylab = 'Dichte', lwd = 2,
      main = 'Parametrische vs. nichtparametrische Dichteschätzung')
lines(density(CRP[!is.na(CRP)], from = 0), lwd = 2)

#######################################
## Nach der Installation des Pakets 'ROptEst' findet sich im Verzeichnis
## .../R/library/ROptEst ein Verzeichnis 'scripts'. Dort findet man weitere
## Beispiele für die Berechnung von AL-Schätzern.
#######################################