Rmd

Lösungen mit angegeben

Die Variable multi.plot steuert, ob beim Optimieren jeder Zwischenschritt visualisiert wird. Dabei werden sehr viele sehr ähnliche Grafiken produziert, die das Verständnis erleichtern sollen, den Output aber lang werden lassen.

multi.plot = 0 # a constant to switch, if 1 every step of optimization will be plotted

Statistische Modelle und Modelltests

cf R for Data Science, S. 345

The goal of a model is to provide a simple low-dimensional summary of a dataset. In the context of this book we’re going to use models to partition data into patterns and residuals. Strong patterns will hide subtler trends, so we’ll use models to help peel back layers of structure as we explore a dataset.
However, before we can start using models on interesting, real datasets, you need to understand the basics of how models work. For that reason, this chapter of the book is unique because it uses only simulated datasets. These datasets are very simple, and not at all interesting, but they will help you understand the essence of modeling before you apply the same techniques to real data in the next chapter.
There are two parts to a model:
1. First, you define a family of models that express a precise, but generic, pattern that you want to capture. For example, the pattern might be a straight line, or a quadatric curve. You will express the model family as an equation like y = a_1 * x + a_2ory = a_1 * x ^ a_2. Here, x and y are known variables from your data, and a_1 and a_2 are parameters that can vary to capture different patterns.
2. Next, you generate a  tted model by finding the model from the family that is the closest to your data. This takes the generic model family and makes it specific, like y = 3 * x + 7 or y = 9 * x ^ 2.
It’s important to understand that a fitted model is just the closest model from a family of models. That implies that you have the “best” model (according to some criteria); it doesn’t imply that you have a good model and it certainly doesn’t imply that the model is “true.” George Box puts this well in his famous aphorism:
All models are wrong, but some are useful.
It’s worth reading the fuller context of the quote:
Now it would be very remarkable if any system existing in the real world could be exactly represented by any simple model. However, cunningly chosen parsimonious models often do provide remarkably useful approximations. For example, the law PV = RT relating pressure P, volume V and temperature T of an “ideal” gas via a con‐ stant R is not exactly true for any real gas, but it frequently provides a useful approximation and furthermore its structure is informative since it springs from a physical view of the behavior of gas molecules.
For such a model there is no need to ask the question “Is the model true?” If “truth” is to be the “whole truth” the answer must be “No.” The only question of interest is “Is the model illuminating and useful?”
The goal of a model is not to uncover truth, but to discover a simple approximation that is still useful.

Wir können die Parameter eines Modells an Daten anpassen, indem wir ein Modell angeben und dessen Parameter an einen Datensatz anpassen lassen. MLE (Maximum Likelihood Estimation) ist ein Optimierungsproblem. In R gibt es hierzu den Befehl optim(). Mit mle() aus der library(stats4) können wir MLE der Parameter erhalten. Ebenso mit library(bbmle), was nicht so sensibel auf nicht gut passende Startwerte reagiert. Außerdem bieten verschiedene Pakete ebenfalls in ihrem Kontext ML-Schätzungen an. Hier zeigen wir das am Beispiel von Multilevel-Analysen mit library(nlme) und dem Befehl lme().

Probability P ist die Wahrscheinlichkeit von Daten unter der Bedingung bestimmter Parameter p. P(X | p)

Likelihood: Wahrscheinlichkeit von Paramtern unter der Voraussetzung bestimmter Daten L(p | X)

Das Ziel von MLE (maximum likelihood estimation) ist es, die Parameterausprägung(-en) zu finden, die die beobachteten Daten am wahrscheinlichlichsten machen. In MLE sind also die Daten gegeben und die Ausprägungen der Parameter werden gesucht/angepasst (Daten-Analyse).

Probability Knowing parameters -> Prediction of outcome

Likelihood Observation of data -> Estimation of parameters

“Klassisches”, deterministisches Vorgehen: Ordinary Least Squares Method (OLS). Visaulisierung hiervon am Beispiel einer Geraden durch eine Punktwolke:

http://students.brown.edu/seeing-theory/regression/index.html#first

Wahrscheinlichkeit (Probability) und Likelihood bei binomialverteilten Daten

In enger Anlehung an Quelle

Beim Münzwurf bzw. bei vielen Münzwürfen hat das Modell einen Parameter, die Wahrscheinlichkeit, eine Seite, z. B. Kopf zu treffen. p(Kopf) = 0.5. Probabilities addieren sich zu 1. (20% Regenwahrscheinlichkeit => 80% Wahrscheinlichkeit für keinen Regen). Bei unabhängigen Ereignissen ist die Gesamtwahrscheinlichkeit (Verbundwahrscheinlichkeit) gleich dem Produkt der Einzelwahrscheinlichkeiten. Wahrscheinlichkeit, bei zwei Würfen zweimal Kopf zu bekommen ist 0.5 * 0.5 = 0.25.

Beispiel: Bei 100 Münzwürfen sind 57 mal Kopf gefallen. Ist die Münze “gerecht”? n=100 k=57

Dichotomes Ereignis. Verteilung: Binomialverteilung Wahrscheinlichkeitsverteilung:

\[ \frac{n!}{k!(n-k)!}p^k(1-p)^{n-k} \] n = Anzahl der Münzwürfe k = Anzahl der “Treffer” (Kopf) p = Wahrscheinlichkeit, bei jedem einzelnen Münzwurf einen Kopf zu bekommen

Der zweite Teil der Formel hat zu tun mit der Verbundwahrscheinlichkeit der Ereignisse. Voraussetzung ist die Unabhängigkeit der Ereignisse und die Konstanz der Auftretenswahrscheinlichkeit. Unter diesen Voraussetzungen ist die Verbundwahrscheinlichkeit gleich dem Produkt der Einzelwahrscheinlichkeiten. Für den Fall von 9 Würfen mit 4 Treffern ergibt das also: \(p^4(1-p)^5 = {0.5}^4(0.5)^5\) Der erste Teil der Formel hat zu tun mit der Anzahl der möglichen Reihenfolgen, mit denen ein solches Ergebnis zustande kommen kann. Vorstellbare Reihenfolgen können sein:

K, Z, K, K, Z, Z, K, Z, Z.

oder

Z, K, K, Z, K, Z, Z, K, Z.

oder sogar

K, K, K, K, Z, Z, Z, Z, Z.

Jede Permutation hat die gleiche Auftretenswahrscheinlichkeit. \(\frac{n!}{k!(n-k)!}\) ist also für n=9 und k=4

\[ \frac{n!}{k!(n-k)!}p^k(1-p)^{n-k} = \frac{9!}{4!(5-4)!}0.5^4(1-0.5)^{9-4} = 0.246 \]

Probability Knowing parameters -> Prediction of outcome

Likelihood Observation of data -> Estimation of parameters

Likelihood der Daten unter der Annahme, die Wahrscheinlichkeit p sei 0.5

n=100 k=57

\[ L(p=0.5 | data) = \frac{100!}{57!(100-57)!}0.5^{57} * 0.5^{43} = 0.03006864 \]

Likelihood surface

# we define a function that calculates likelihood of k hits in n events with probability of p 
l.f = function(p,n,k){return(factorial(n)/(factorial(k) * factorial(n-k)) * p^k * (1-p)^(n-k))}
# L with p=0.5, n=100, k=57
l.f(0.5, 100,57)
## [1] 0.03006864
# L with p=0.52, n=100, k=56
l.f(0.52,100,57)
## [1] 0.04860458
# we generate a vector of probabilities
probs <- seq(0.01, 0.99, by=0.01)
# and calculate the likelihoods for these probabilities and n=100 and k=57
lik.probs <- l.f(probs, 100, 57)
plot(probs, lik.probs)

# consider: a=b*c then log(a) = log(b) + log(c)
plot(probs, log(lik.probs))

plot(probs, -2*log(lik.probs))

Die Multiplikation vieler kleiner Probabilities erzeugt sehr sehr kleine Resultate. Logarithmieren hilft hier. Der Logarithmus von Werten zwischen 0 und 1 ist immer negativ. Die Multiplikation von Probabilities ändert sich zur Addition(!) ihrer Logarithmen. \(a=b*c =>log(a) = log(b) + log(c)\)

a=12;b=4; c=3
b*c
## [1] 12
log(a)
## [1] 2.484907
log(b) + log(c)
## [1] 2.484907
# exp() ist Umkehrfunktion von log()
exp(log(a))
## [1] 12
exp(log(b) + log(c))
## [1] 12
# we define a function that calculates likelihood of k hits in n events with probability of p 
l.f = function(p,n,k){return(factorial(n)/(factorial(k) * factorial(n-k)) * p^k * (1-p)^(n-k))}

# let's calculate 5 events: 3 times head and 2 times tail
# probablity
p.5 = factorial(5)/(factorial(3) * factorial(5-3)) * 0.5^3 * (1 - 0.5)^(5 - 3)  
p.5
## [1] 0.3125
l.5 <- log(p.5)
l.5
## [1] -1.163151

likelihood ratio test beim Münzwurf

Beim Likelihood Ratio Test werden zwei Likelihoods zueinander in Beziehung gesetzt. Dies sind normalerweise geschachtelte Modelle, wobei das komplexere Modell gegen ein einfacheres Modell mit weniger Parametern, als das komplexere Modell getestet wird.

In unserem Münz-Beispiel:

Alternativhypothese H1: p ist ungleich 0.5

Nullhypothese H0: p = 0.5

Der Likelihood Ratio Test beantwortet die Frage: Sind die Daten signifikant weniger wahrscheinlich unter der Nullhypothese entstanden, als unter der Alternativhypothese.

Anders ausgedrückt: Fitted das komplexere Modell mit mehr Parametern die Daten signifikant besser, als das einfachere Modell? “Lohnt” sich die Aufnahme der zusätzlichen Parameter im komplexeren Modell?

Vorgehen: Berechnen der Likelihood Ratio:

Differenz der Log-Likelihoods zwischen den beiden Modellen. Da a=b/c => log(a) = log(b) - log(c) sprechen wir von einer Ratio (Bruch), obwohl wir de facto eine Differenz bilden.

Prüfgröße ( \(\chi^2\) verteilt): \(2 ( LL_A - LL_0)\)

Die Multiplikation mit 2 erfolgt, weil die Differenz der Log-Likelihoods \(\chi^2\)-verteilt ist.

Die Freiheitsgrade für den Test ergeben sich aus der Differenz der Parameter im Alternativ- zum Nullmodell. In unserem Beispiel hat die Alternativhypothese 1 Parameter (p), der geschätzt wird und das Nullmodell hat keinen Parameter (p ist gegeben).

Maß H1 H0
p 0.57 0.50
Likelihood 0.08037551 0.0389
Log Likelihood -2.521046 -3.247

\[ 2(LL_A - LL_0) = 2 * (-2.521046 + 3.247) = 1.451908 \]

Kritischer \(\chi^2\)-Wert bei einem Freiheitsgrad ist 3.84 (qchisq(0.95, df=1) -> 3.8414588) Der empirische \(\chi^2\)-Wert ist also kleiner als der kritische \(\chi^2\)-Wert. Wir schließen, dass der Model-Fit unter der Alternativhypothese nicht sig. schlechter ist, als unter der Nullhypothese. Wir haben also keinen Grund, an der “Gerechtigkeit” der Münze zu zweifeln.

“klassischer” Binomialtest

Nur zur Ergänzung …

binom.test(57, 100, p = 0.5)
## 
##  Exact binomial test
## 
## data:  57 and 100
## number of successes = 57, number of trials = 100, p-value = 0.1933
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.4671337 0.6686090
## sample estimates:
## probability of success 
##                   0.57

Aufgabe zum Münzwurf

  • Berechnen Sie analog zum obigen Beispiel die Likelihood von 84 mal Zahl geworfen zu haben bei 123 Durchgängen, wenn Zahl das gewünschte Ereignis ist.
  • Berechnen Sie die -2LL des Experiments
  • Prüfen Sie mit dem Likelihood-Ratio-Test, ob die Münze in diesem Experiment “gerecht” geblieben ist.

Lösung

# n is 123, k is 84
# we define a function that calculates likelihood of k hits in n events with probability of p 
l.f = function(p,n,k){return(factorial(n)/(factorial(k) * factorial(n-k)) * p^k * (1-p)^(n-k))}
probs <- seq(0.01, 0.99, by=0.01)
lik.probs <- l.f(probs, 123, 84)
max.lik <- max(lik.probs)
max.lik
## [1] 0.07692596
# let's calculate 123 events: 84 times head and 39 times tail
# probablity
pp = factorial(84)/(factorial(84) * factorial(123 - 84)) * 0.5^84 * (1 - 0.5)^(123 - 84)  
pp
## [1] 4.61026e-84
pp = l.f(max.lik, 123, 84)
pp
## [1] 2.132192e-63
p.50 = l.f(.5, 123, 84)
emp.chisquare <- -2 * log(pp) - (-2 * (log(p.50)))
emp.chisquare
## [1] 266.6346
crit.chisquare <- qchisq(0.95, df=1)
crit.chisquare
## [1] 3.841459
if (emp.chisquare > crit.chisquare) {cat("sig")} else {cat("n.s.")}
## sig

Würfelbeispiel

Probability bei Würfeln Likelihood beim Würfeln

Bei einem Würfel wird vermutet, er sei manipuliert. Frage: Ist ein Würfel manipuliert? 132 mal gewürfelt, 36 mal fällt die 6.

\[ L(\pi|x_1, x_2, ..., x_n) = {132 \choose 36} * \pi^{36} * (1-\pi)^{96} \] wobei

\[ {n \choose k} = \frac{n!}{k! * (n-k)!} \]

\[ \frac{dL}{d\pi} = {132 \choose 36}36\pi^{35}(1-\pi)^{96}-{132 \choose 36}96\pi^{36}(1-\pi)^{95} = 0 \] \[ => 3(1-\pi)-8\pi = 0 \] \[ => 3 - 3\pi - 8\pi = 0 \] \[ => 11 \pi = 3 \] \[ => \hat{\pi} = 3/11 = 0.27... \]

Aufgaben

  • Berechnen Sie die Likelihood von 30 5en bei 120 Würfelwürfen.
  • Berechnen Sie die Likelihood von 45 3en bei 180 Würfelwürfen.
  • Was fällt auf?

Lösung

# l.f works for a maximum of n=170 only because factorial() works only to this limit
l.f = function(p,n,k){if (n < 171) {return(factorial(n)/(factorial(k) * factorial(n-k)) * p^k * (1-p)^(n-k))} else {return(NA)}}
lik.ratio.test = function(p, n, k){
  probs <- seq(0.01, 0.99, by=0.01)
  lik.probs <- l.f(probs, n, k)
  max.lik <- max(lik.probs)
  cat('max.lik: ')
  cat(max.lik)
  # let's calculate 123 events: 84 times head and 39 times tail
  # probablity
  pp = factorial(k)/(factorial(k) * factorial(n - k)) * p^k * (1 - p)^(n-k)  
  # pp
  pp = l.f(max.lik, n, k)
  cat(' data lik: ')
  cat(pp)
  pc = l.f(p, n, k)
  emp.chisquare <- -2 * log(pp) - (-2 * (log(pc)))
  #emp.chisquare
  crit.chisquare <- qchisq(0.95, df=1)
  #crit.chisquare
  if (emp.chisquare > crit.chisquare) {cat(" sig ")} else {cat(" n.s. ")}
  return(paste('emp.chisquare: ', emp.chisquare, ' crit.chisquare: ', crit.chisquare))
}

lik.ratio.test(p=1/6, n=120, k=30)
## max.lik: 0.08385171 data lik: 3.251442e-08 sig
## [1] "emp.chisquare:  24.1627400338251  crit.chisquare:  3.84145882069412"
# this would throw an error because n > 170, so it is commented out
# lik.ratio.test(p=1/6, n=180, k=43)


# control
binom.test(30, 120, p = 1/6)
## 
##  Exact binomial test
## 
## data:  30 and 120
## number of successes = 30, number of trials = 120, p-value = 0.01933
## alternative hypothesis: true probability of success is not equal to 0.1666667
## 95 percent confidence interval:
##  0.1754646 0.3372692
## sample estimates:
## probability of success 
##                   0.25
binom.test(43, 180, p = 1/6)
## 
##  Exact binomial test
## 
## data:  43 and 180
## number of successes = 43, number of trials = 180, p-value = 0.01217
## alternative hypothesis: true probability of success is not equal to 0.1666667
## 95 percent confidence interval:
##  0.1786013 0.3079847
## sample estimates:
## probability of success 
##              0.2388889

Likelihood bei normalverteilten Daten

Eine Normalverteilung wird vollständig beschrieben durch Mittelwert und Standardabweichung, also durch die beiden Parameter \(\bar{x}, sd\).

Ermittlung der Likelihood für gegebene Daten

Die Likelihood normalverteilter Daten ist gleich der Multiplikation aller Einzel-Wahrscheinlichkeiten. Die Einzel-Wahrscheinlichkeiten können wir in R über die Funktion dnorm() erhalten.

rnorm(100, 0, 2) erzeugt 100 normalverteilte Zufallswerte \(x_1, ..., x_100\) mit \(\bar{x}\) und \(sd=2\)

dnorm(1, 0, 2) Probability Density Function (PDF) Höhe der Wahrscheinlichkeitsfunktion an einer bestimmten Stelle, hier 1, für eine Normalverteilung mit \(\bar{x}\) und \(sd=2\)

pnorm(1.96,0,2) Cumulative Distribution Function (CDF) gibt die Fläche unter der NV von \(-\infty\) bis 1.96

qnorm(0.975,0,2) Quantile Function – Inverse von pnorm(), gibt den x-Wert, an dem die CDF (cumulative Fläche) der NV(0,2) den Wert 0.975 annimmt.

# we take a first look at standard normal distribution in this context
# we guarantee reproducable random values
set.seed(42)
dd <- rnorm(100, 0, 1)
head(cbind(dd, dnorm(dd, 0, 1)))
##              dd          
## [1,]  1.3709584 0.1558748
## [2,] -0.5646982 0.3401459
## [3,]  0.3631284 0.3734879
## [4,]  0.6328626 0.3265422
## [5,]  0.4042683 0.3676386
## [6,] -0.1061245 0.3967021
# first two probabilities
dnorm(dd, 0, 1)[1]
## [1] 0.1558748
dnorm(dd, 0, 1)[2]
## [1] 0.3401459
# likelihood is the product of all single probabilities
# first two elements
dnorm(dd, 0, 1)[1] * dnorm(dd, 0, 1)[2]
## [1] 0.05302018
# multiplication of probabilities results in really small values
prod(dnorm(dd, 0, 1))
## [1] 5.695794e-64
# logarithm of it makes it more handy
log(dnorm(dd, 0, 1)[1] * dnorm(dd, 0, 1)[2])
## [1] -2.937083
# using log of probabilities and sum it up is equivalent
log(dnorm(dd, 0, 1)[1])
## [1] -1.858702
log(dnorm(dd, 0, 1)[2])
## [1] -1.078381
log(dnorm(dd, 0, 1)[1]) + log(dnorm(dd, 0, 1)[2])
## [1] -2.937083
# inverse should result in the same number
exp(log(dnorm(dd, 0, 1)[1] * dnorm(dd, 0, 1)[2]))
## [1] 0.05302018
exp(log(dnorm(dd, 0, 1)[1]) + log(dnorm(dd, 0, 1)[2]))
## [1] 0.05302018
# so we end up calculating -log(likelihood) using
-sum(log(dnorm(dd, 0, 1)))
## [1] 145.6257
# or the often used form of -2LL
-2 * sum(log(dnorm(dd, 0, 1)))
## [1] 291.2514

Likelihood-Quotienten-Test (Likelihood Ratio Test) bei normalverteilten Daten

Wenn wir einen Datensatz angeblich normalverteilter Werte haben, können wir nun die Likelihood berechnen, die wir unter der Annahme bestimmter NV-Parameter erhalten würden.

Angenommen, wir wüssten vom unten berechneten Datenvektor nur, dass er normalverteilt sei. Wir vermuten, die Daten kommen aus einer Standard-Normalverteilung N(0,1). Dies ist unser Nullmodell. Ein Kollege “schätzt” den Mittelwert aus den gegebenen Daten auf etwa 1 und “verbraucht” hierfür einen Freiheitsgrad. Wir wollen zunächst wissen, ob der nun geschätzte Mittelwert die Daten besser erklärt, als das Nullmodell, also, ob der Fit besser ist. Als Test benutzen wir den Likelihood-Quotienten-Test

\[ 2(LL_{Alternativmodell} - LL_{Nullmodell})) \]

Die resultierende Prüfgröße ist \(\chi^2\)-verteilt

Unser Nullmodell hat 0 Freiheitsgrade, das Alternativmodell 1 df. Kritischer \(\chi^2\)-Wert bei einem Freiheitsgrad ist 3.84 (qchisq(0.95, df=1) -> 3.8414588)

# we guarantee reproducable random values
set.seed(42)
dd <- rnorm(100, 1.2, 2)
# we add some noise
set.seed(43)
dd <- dd + rnorm(100, 0, 0.1)

# let's assume, we only have data dd and the information, that it's normally distributed data, nothing more
# likelihood that these data come from a standard normal distribution N(0,1)
LL.0.1 <- sum(log(dnorm(dd, 0, 1)))
LL.0.1
## [1] -387.4276
# some colleague might have informed us, that he "estimated from the data", that the data mean should be around 1 and he used one df for that
# so we calculate likelihood of N(1,1)
LL.1.1 <- sum(log(dnorm(dd, 1, 1)))
LL.1.1
## [1] -310.3015
# clearly, log likelihood of the second model fits the data better (is less negative)
# but is this significant?
# we use likelihood ratio test
# chi^2 test statistik
2 * (LL.1.1 - LL.0.1)
## [1] 154.2522
# critical chi^2 value for 1 df
qchisq(0.95, df=1)
## [1] 3.841459
# test statistik > critical value?
2 * (LL.1.1 - LL.0.1) > qchisq(0.95, df=1)
## [1] TRUE
# empirical value is bigger than critical value
# so our alternative model fits the data significantly better than the zero model

Aber: Alle Infos zu obiger Prüfung kamen “von außen”

Nun prüfen wir, ob eine MLE (maximum likelihood estimation) des Mittelwertes und der Standardabweichung das Modell sig. verbessert

Idee: Wir vergleichen ein Modell, in dem wir MLE Schätzungen von MW und sd vergleichen mit dem Modell des Kollegen, wo nur der Mittelwert “geschätzt” worden war.

Vorgehen: Wir brauchen ein Verfahren, mit dem wir die Parameter einer NV aus vorgegebenen Daten schätzen können.

Anpassen einer Normalverteilung: MLE des MW und sd einer NV

cf

Wir benutzen die Funktion mle() des Pakets stats4. Um Parameter schätzen zu können, braucht mle() eine Funktion und sowie die Parameter, die geschätzt werden sollen.

# we set random seed to get reproducalble random data
set.seed(42)
N <- 100
x <- rnorm(N, mean = 3, sd = 2)
mean(x)
## [1] 3.06503
sd(x)
## [1] 2.082714
# Then we formulate the log-likelihood function.
LL <- function(mu, sigma) {
    R = dnorm(x, mu, sigma)
    -sum(log(R))
}

# And apply MLE to estimate the two parameters (mean and standard deviation) for which the normal distribution best describes the data.
require(stats4)
## Loading required package: stats4
mle(LL, start = list(mu = 1, sigma=1))
## Warning in dnorm(x, mu, sigma): NaNs wurden erzeugt

## Warning in dnorm(x, mu, sigma): NaNs wurden erzeugt

## Warning in dnorm(x, mu, sigma): NaNs wurden erzeugt

## Warning in dnorm(x, mu, sigma): NaNs wurden erzeugt

## Warning in dnorm(x, mu, sigma): NaNs wurden erzeugt

## Warning in dnorm(x, mu, sigma): NaNs wurden erzeugt

## Warning in dnorm(x, mu, sigma): NaNs wurden erzeugt
## 
## Call:
## mle(minuslogl = LL, start = list(mu = 1, sigma = 1))
## 
## Coefficients:
##       mu    sigma 
## 3.065028 2.072280
# We don't want to see the warnings, produced when negative values are attempted for the standard deviation.
LL <- function(mu, sigma) {
    R = suppressWarnings(dnorm(x, mu, sigma))
    #
    -sum(log(R))
    }
mle(LL, start = list(mu = 1, sigma=1))
## 
## Call:
## mle(minuslogl = LL, start = list(mu = 1, sigma = 1))
## 
## Coefficients:
##       mu    sigma 
## 3.065028 2.072280

Zurück zu obigem Beispiel und der Frage, ob es eine besser angepasste NV für unsere Daten gibt, als die, die der Kollege vorschlug: N(1,1)

# we guarantee reproducable random values
set.seed(42)
dd <- rnorm(100, 1.2, 2)
# we add some noise
set.seed(43)
dd <- dd + rnorm(100, 0, 0.1)

# LL of the colleague's model N(1,1)
LL.1.1 <- sum(log(dnorm(dd, 1, 1)))
LL.1.1
## [1] -310.3015
# we make MLE for mean and sd
# we define an adequate function to optimize
LL <- function(mu, sigma) {
    R = suppressWarnings(dnorm(dd, mu, sigma))
    return(-sum(log(R)))
    }

# now we do MLE
p.mle <- mle(LL, start = list(mu = 1, sigma=1))
# mle of our parameters
p.mle@coef
##       mu    sigma 
## 1.271261 2.072335
# we now calculate the likelihood of the data for our fitted normal distribution
LL.mle <- sum(log(dnorm(dd, p.mle@coef[1], p.mle@coef[2])))
LL.mle
## [1] -214.7614
# chi^2 test statistik
2 * (LL.mle - LL.1.1)
## [1] 191.0802
# critical chi^2 value for 1 df
qchisq(0.95, df=1)
## [1] 3.841459
# test statistik > critical value?
2 * (LL.mle - LL.1.1) > qchisq(0.95, df=1)
## [1] TRUE
# empirical value is bigger than critical value
# so our alternative model fits the data significantly better than the zero model

Übung:

  • generieren Sie eine NV mit 100 Werten, einem Mittelwert von 108 und einer Standardabweichung von 12.5
  • dies seien Daten, von denen Sie nur wissen, dass sie aus einer Intelligenzmessung stammen und in der IQ-Norm vorliegen: N(100,15)
  • ermitteln Sie die Likelihood der vorliegenden Daten unter der Annahme, es seien IQs
  • führen Sie eine MLE der Paramter der vorliegenden Daten durch unter der Annahme, die Daten seien normalverteilt
  • prüfen Sie, ob die von Ihnen ermittelten Parameter die Daten besser erklären, als die Annahme, es seien IQ-Daten
  • welche Schlüsse ziehen Sie aus Ihren Berechnungen
  • überlegen Sie sich ein alternatives (“klassisches”) Vorgehen zur Überprüfung obiger Fragestellung

Verstehbeispiel einfache Regression mit Parameterschätzung über verschiedene Methoden

  • Suche nach den Parametern nach der least square Methode
  • Suche nach den Parameter nach der maximum likelihood Methode.
  • Vergleich mit den Parametern eines “klassischen” lm().

Wir brauchen ein generisches Modell (model family). Hier: g.model() Es berechnet nach Parameterübergabe die geschätzten Werte für die Reaktionsvariable.

Für die Optimierung mit dem Befehl optim() brauchen wir eine Funktion, deren erster Parameter ein Array mit Werten für die Modellparameter sind. Hier sind das die beiden Funktionen measure_distance() für die Methode der Quadratsummen-Minimierung und m.ml() für die Maximum-Likelihood-Schätzung. Bei der Quadratsummen-Minimierung werden die von optim() übergebenen Parameter so lange optimiert, bis der Rückgabewert, die Abweichungsquadrate der Schätzwerte von den Werten der Reaktionsvariable minimal sind. Bei der ML-Schätzung werden die von optim() übergebenen Parameter so lange optimiert, bis der zurückgegebene Likelihood-Wert minimal wird. Für die Ermittlung der Likelihood wird die Density der Residuen berechnet, die sich um einen Mittelwert von 0 normal verteilen (müssen). Als sd der NV wird hierbei die sd der Differenzen y_geschätzt - y eingesetzt.

Least Square visualisiert

Selbst spielen mit einer Shiny-App von Andreas Cordes

# some data
df <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/ml_lm_simple_1.txt")
# we define a function that implements the model family, a generic model
g.model <- function(a, data) {
  a[1] + data$x * a[2]
    }
ii <- 0
head(g.model(c(0, 2), df))
## [1]  4.591114  6.370945 -1.585213 -1.900441 -2.866911  2.154753
require(ggplot2)
## Loading required package: ggplot2
# another function to calculate the square mean distances which is the criterion to optimize
measure_distance <- function(mod, data) {
    diff <- data$y - g.model(mod, data)
    if (multi.plot) {
        ii <<- ii + 1
        pp <- ggplot(data=data, aes(x=x, y=y)) + geom_point() + geom_abline(intercept = mod[1], slope = mod[2], color="red")
        ggsave(paste("lsq_", formatC(ii, width=3, flag="0"), ".png", sep=""), units="in", width=8, height=6, dpi=72)
        }
    print(paste(ii, round(mod[1], 3), round(mod[2],3), round(sqrt(mean(diff ^2)), 3)))
    return(sqrt(mean(diff ^ 2)))
}
measure_distance(c(0,2), df) 
## [1] "0 0 2 10.122"
## [1] 10.12201
ii <- 0
best.lsq <- optim(c(0, 0), measure_distance, data = df) 
## [1] "0 0 0 10.535"
## [1] "0 0.1 0 10.444"
## [1] "0 0 0.1 10.468"
## [1] "0 0.1 0.1 10.376"
## [1] "0 0.15 0.15 10.297"
## [1] "0 0.25 0.05 10.272"
## [1] "0 0.375 0.025 10.175"
## [1] "0 0.425 0.175 10.027"
## [1] "0 0.588 0.263 9.821"
## [1] "0 0.813 0.138 9.696"
## [1] "0 1.144 0.131 9.399"
## [1] "0 1.356 0.369 9.045"
## [1] "0 1.847 0.541 8.49"
## [1] "0 2.403 0.409 8.062"
## [1] "0 3.311 0.483 7.195"
## [1] "0 4.014 0.892 6.306"
## [1] "0 5.449 1.273 4.849"
## [1] "0 6.913 1.215 3.617"
## [1] "0 9.446 1.552 1.995"
## [1] "0 11.585 2.342 3"
## [1] "0 9.516 1.877 2.048"
## [1] "0 13.514 2.156 4.249"
## [1] "0 11.497 1.935 2.606"
## [1] "0 7.465 1.494 3.086"
## [1] "0 10.489 1.825 2.078"
## [1] "0 8.473 1.604 2.384"
## [1] "0 9.985 1.77 1.967"
## [1] "0 9.916 1.445 1.998"
## [1] "0 9.816 1.553 1.955"
## [1] "0 10.355 1.771 2.023"
## [1] "0 9.674 1.607 1.953"
## [1] "0 9.504 1.39 2.053"
## [1] "0 9.865 1.675 1.944"
## [1] "0 9.723 1.729 1.958"
## [1] "0 9.793 1.597 1.947"
## [1] "0 9.984 1.665 1.948"
## [1] "0 9.906 1.65 1.944"
## [1] "0 9.979 1.728 1.956"
## [1] "0 9.839 1.63 1.943"
## [1] "0 9.88 1.605 1.945"
## [1] "0 9.869 1.657 1.943"
## [1] "0 9.802 1.637 1.944"
## [1] "0 9.88 1.647 1.943"
## [1] "0 9.91 1.675 1.945"
## [1] "0 9.857 1.641 1.943"
## [1] "0 9.868 1.63 1.943"
## [1] "0 9.869 1.651 1.943"
## [1] "0 9.845 1.645 1.943"
## [1] "0 9.854 1.645 1.943"
## [1] "0 9.842 1.635 1.943"
## [1] "0 9.862 1.647 1.943"
## [1] "0 9.859 1.651 1.943"
## [1] "0 9.857 1.643 1.943"
## [1] "0 9.849 1.642 1.943"
## [1] "0 9.859 1.646 1.943"
## [1] "0 9.853 1.643 1.943"
## [1] "0 9.857 1.645 1.943"
## [1] "0 9.854 1.647 1.943"
## [1] "0 9.857 1.644 1.943"
## [1] "0 9.86 1.644 1.943"
## [1] "0 9.856 1.645 1.943"
## [1] "0 9.855 1.644 1.943"
## [1] "0 9.857 1.645 1.943"
best.lsq$par
## [1] 9.856692 1.644704
plot(df$x, df$y)
abline(coef=best.lsq$par, col = "red")

# we now estimate constant and slope using maximum likelihood estimation
m.ml <- function(mod, data) {
        yHat = g.model(mod, data)
        error = data$y - yHat
        output = -sum(log(dnorm(error, mean = 0, sd=sd(error))))
        if (multi.plot) {
      ii <<- ii + 1
      pp <- ggplot(data=data, aes(x=x, y=y)) + geom_point() + geom_abline(intercept = mod[1], slope = mod[2], color="red")
            ggsave(paste("ml_", formatC(ii, width=3, flag="0"), ".png", sep=""), units="in", width=8, height=6, dpi=72)
      #ggsave(paste("ml_", ii, ".png", sep=""))
            # plot(data$x,data$y)
            # abline(coef=c(mod[1], mod[2]), col = "red")
            }
        print(paste(ii, "-ll: ", output, " b0: ", round(mod[1], 3), " b1: ", round(mod[2], 3), " mean err: ", round(mean(error), 3), " sd.err: ", round(sd(error), 3)))
        return(output)
        }

m.ml(c(0,1), df)
## [1] "0 -ll:  1027.37595222304  b0:  0  b1:  1  mean err:  9.776  sd.err:  2.449"
## [1] 1027.376
ii <- 0
best.ml <- optim(c(0, 1), m.ml, data = df) 
## [1] "0 -ll:  1027.37595222304  b0:  0  b1:  1  mean err:  9.776  sd.err:  2.449"
## [1] "0 -ll:  1011.16562531869  b0:  0.1  b1:  1  mean err:  9.676  sd.err:  2.449"
## [1] "0 -ll:  1116.8666581843  b0:  0  b1:  1.1  mean err:  9.788  sd.err:  2.318"
## [1] "0 -ll:  930.392774078897  b0:  0.1  b1:  0.9  mean err:  9.663  sd.err:  2.594"
## [1] "0 -ll:  851.63967368255  b0:  0.15  b1:  0.8  mean err:  9.601  sd.err:  2.751"
## [1] "0 -ll:  839.01774908585  b0:  0.25  b1:  0.8  mean err:  9.501  sd.err:  2.751"
## [1] "0 -ll:  763.644891986868  b0:  0.375  b1:  0.7  mean err:  9.363  sd.err:  2.917"
## [1] "0 -ll:  662.813770230545  b0:  0.425  b1:  0.5  mean err:  9.288  sd.err:  3.272"
## [1] "0 -ll:  567.914755495011  b0:  0.588  b1:  0.25  mean err:  9.094  sd.err:  3.748"
## [1] "0 -ll:  530.608649650061  b0:  0.813  b1:  0.15  mean err:  8.857  sd.err:  3.945"
## [1] "0 -ll:  463.699078099855  b0:  1.144  b1:  -0.175  mean err:  8.485  sd.err:  4.608"
## [1] "0 -ll:  422.132444819916  b0:  1.356  b1:  -0.625  mean err:  8.216  sd.err:  5.56"
## [1] "0 -ll:  395.577811755424  b0:  1.847  b1:  -1.288  mean err:  7.643  sd.err:  7.003"
## [1] "0 -ll:  387.828938603691  b0:  2.403  b1:  -1.713  mean err:  7.033  sd.err:  7.944"
## [1] "0 -ll:  390.572177938474  b0:  3.311  b1:  -2.694  mean err:  6.003  sd.err:  10.141"
## [1] "0 -ll:  393.520847442485  b0:  3.106  b1:  -2.825  mean err:  6.191  sd.err:  10.437"
## [1] "0 -ll:  389.113924633951  b0:  2.616  b1:  -2.163  mean err:  6.764  sd.err:  8.948"
## [1] "0 -ll:  389.987328121072  b0:  3.172  b1:  -2.588  mean err:  6.155  sd.err:  9.902"
## [1] "0 -ll:  388.331709939349  b0:  2.841  b1:  -2.263  mean err:  6.527  sd.err:  9.172"
## [1] "0 -ll:  386.022728842732  b0:  2.628  b1:  -1.813  mean err:  6.796  sd.err:  8.167"
## [1] "0 -ll:  384.867337708838  b0:  2.634  b1:  -1.638  mean err:  6.811  sd.err:  7.778"
## [1] "0 -ll:  391.694510669089  b0:  2.197  b1:  -1.088  mean err:  7.318  sd.err:  6.564"
## [1] "0 -ll:  386.762490022437  b0:  2.68  b1:  -1.969  mean err:  6.725  sd.err:  8.515"
## [1] "0 -ll:  383.934662759193  b0:  2.911  b1:  -1.894  mean err:  6.503  sd.err:  8.348"
## [1] "0 -ll:  382.598921611552  b0:  3.165  b1:  -1.984  mean err:  6.238  sd.err:  8.55"
## [1] "0 -ll:  379.730649661993  b0:  3.12  b1:  -1.653  mean err:  6.324  sd.err:  7.812"
## [1] "0 -ll:  376.055472363282  b0:  3.339  b1:  -1.495  mean err:  6.124  sd.err:  7.462"
## [1] "0 -ll:  374.932095404623  b0:  3.87  b1:  -1.842  mean err:  5.55  sd.err:  8.233"
## [1] "0 -ll:  371.84630993749  b0:  4.488  b1:  -1.945  mean err:  4.92  sd.err:  8.461"
## [1] "0 -ll:  362.432109674035  b0:  4.662  b1:  -1.455  mean err:  4.806  sd.err:  7.374"
## [1] "0 -ll:  351.096157938393  b0:  5.411  b1:  -1.191  mean err:  4.091  sd.err:  6.791"
## [1] "0 -ll:  353.471161651703  b0:  6.559  b1:  -1.64  mean err:  2.886  sd.err:  7.784"
## [1] "0 -ll:  357.066141205344  b0:  5.754  b1:  -1.604  mean err:  3.696  sd.err:  7.703"
## [1] "0 -ll:  328.284190544333  b0:  7.483  b1:  -0.887  mean err:  2.057  sd.err:  6.126"
## [1] "0 -ll:  302.943081171015  b0:  8.98  b1:  -0.358  mean err:  0.626  sd.err:  4.991"
## [1] "0 -ll:  291.725953476172  b0:  7.832  b1:  0.091  mean err:  1.83  sd.err:  4.063"
## [1] "0 -ll:  246.873290101847  b0:  8.468  b1:  0.957  mean err:  1.302  sd.err:  2.51"
## [1] "0 -ll:  269.348198332873  b0:  12.037  b1:  1.79  mean err:  -2.163  sd.err:  1.981"
## [1] "0 -ll:  231.606983340821  b0:  10.381  b1:  1.045  mean err:  -0.599  sd.err:  2.388"
## [1] "0 -ll:  235.08813040406  b0:  9.869  b1:  2.36  mean err:  0.077  sd.err:  2.551"
## [1] "0 -ll:  209.012783778657  b0:  9.647  b1:  1.681  mean err:  0.214  sd.err:  1.955"
## [1] "0 -ll:  245.927401968368  b0:  11.559  b1:  1.768  mean err:  -1.687  sd.err:  1.973"
## [1] "0 -ll:  220.237398532161  b0:  10.786  b1:  1.566  mean err:  -0.94  sd.err:  1.961"
## [1] "0 -ll:  226.266856646731  b0:  10.052  b1:  2.201  mean err:  -0.126  sd.err:  2.333"
## [1] "0 -ll:  213.749092351196  b0:  10.134  b1:  1.912  mean err:  -0.245  sd.err:  2.047"
## [1] "0 -ll:  226.542513189433  b0:  8.995  b1:  2.027  mean err:  0.909  sd.err:  2.141"
## [1] "0 -ll:  211.400543458054  b0:  10.338  b1:  1.681  mean err:  -0.478  sd.err:  1.955"
## [1] "0 -ll:  210.886121083335  b0:  9.851  b1:  1.449  mean err:  -0.019  sd.err:  2.003"
## [1] "0 -ll:  208.829790110259  b0:  9.922  b1:  1.565  mean err:  -0.075  sd.err:  1.961"
## [1] "0 -ll:  213.700048768917  b0:  9.23  b1:  1.565  mean err:  0.616  sd.err:  1.961"
## [1] "0 -ll:  208.873028961433  b0:  10.061  b1:  1.652  mean err:  -0.204  sd.err:  1.953"
## [1] "0 -ll:  212.265262573043  b0:  10.336  b1:  1.536  mean err:  -0.494  sd.err:  1.969"
## [1] "0 -ll:  208.341190335498  b0:  9.819  b1:  1.645  mean err:  0.037  sd.err:  1.953"
## [1] "0 -ll:  209.197503176366  b0:  9.679  b1:  1.558  mean err:  0.166  sd.err:  1.963"
## [1] "0 -ll:  208.504176435714  b0:  9.966  b1:  1.628  mean err:  -0.112  sd.err:  1.953"
## [1] "0 -ll:  208.598672079224  b0:  9.863  b1:  1.708  mean err:  0.001  sd.err:  1.958"
## [1] "0 -ll:  208.379943642962  b0:  9.878  b1:  1.672  mean err:  -0.018  sd.err:  1.954"
## [1] "0 -ll:  208.678434811298  b0:  9.731  b1:  1.688  mean err:  0.131  sd.err:  1.955"
## [1] "0 -ll:  208.357253546875  b0:  9.907  b1:  1.643  mean err:  -0.051  sd.err:  1.953"
## [1] "0 -ll:  208.380244647026  b0:  9.848  b1:  1.616  mean err:  0.004  sd.err:  1.954"
## [1] "0 -ll:  208.337715162232  b0:  9.87  b1:  1.658  mean err:  -0.013  sd.err:  1.953"
## [1] "0 -ll:  208.41312161359  b0:  9.782  b1:  1.659  mean err:  0.076  sd.err:  1.953"
## [1] "0 -ll:  208.328482036305  b0:  9.876  b1:  1.647  mean err:  -0.019  sd.err:  1.953"
## [1] "0 -ll:  208.403892647517  b0:  9.927  b1:  1.661  mean err:  -0.069  sd.err:  1.953"
## [1] "0 -ll:  208.325681176616  b0:  9.846  b1:  1.649  mean err:  0.011  sd.err:  1.953"
## [1] "0 -ll:  208.326244836952  b0:  9.851  b1:  1.638  mean err:  0.004  sd.err:  1.953"
## [1] "0 -ll:  208.323207595383  b0:  9.856  b1:  1.643  mean err:  0  sd.err:  1.953"
## [1] "0 -ll:  208.334619331201  b0:  9.826  b1:  1.644  mean err:  0.03  sd.err:  1.953"
## [1] "0 -ll:  208.323976882781  b0:  9.863  b1:  1.647  mean err:  -0.007  sd.err:  1.953"
## [1] "0 -ll:  208.328187459294  b0:  9.874  b1:  1.641  mean err:  -0.018  sd.err:  1.953"
## [1] "0 -ll:  208.323518035945  b0:  9.853  b1:  1.647  mean err:  0.003  sd.err:  1.953"
## [1] "0 -ll:  208.324570620392  b0:  9.846  b1:  1.643  mean err:  0.01  sd.err:  1.953"
## [1] "0 -ll:  208.323226188336  b0:  9.859  b1:  1.646  mean err:  -0.003  sd.err:  1.953"
## [1] "0 -ll:  208.324043925689  b0:  9.862  b1:  1.642  mean err:  -0.006  sd.err:  1.953"
## [1] "0 -ll:  208.323114066842  b0:  9.855  b1:  1.646  mean err:  0.001  sd.err:  1.953"
## [1] "0 -ll:  208.323406945204  b0:  9.852  b1:  1.643  mean err:  0.003  sd.err:  1.953"
## [1] "0 -ll:  208.323066029986  b0:  9.857  b1:  1.645  mean err:  -0.001  sd.err:  1.953"
## [1] "0 -ll:  208.323652216232  b0:  9.856  b1:  1.648  mean err:  0  sd.err:  1.953"
## [1] "0 -ll:  208.323048337163  b0:  9.856  b1:  1.644  mean err:  0  sd.err:  1.953"
## [1] "0 -ll:  208.323174072351  b0:  9.858  b1:  1.644  mean err:  -0.002  sd.err:  1.953"
## [1] "0 -ll:  208.323051653908  b0:  9.856  b1:  1.645  mean err:  0  sd.err:  1.953"
## [1] "0 -ll:  208.323062514675  b0:  9.855  b1:  1.644  mean err:  0.001  sd.err:  1.953"
## [1] "0 -ll:  208.323041337059  b0:  9.856  b1:  1.644  mean err:  0.001  sd.err:  1.953"
# require(animation)
# #imconvertstring<-'\"/Applications/ImageMagick-7.0.5/bin/convert\" -delay 1x%d %s*.png %s.%s'
# imconvertstring<-'\"/Applications/ImageMagick-7.0.5/bin/convert\"'
# 
# #imconvertstring<-"\"c:\\Program Files\\ImageMagick-6.9.0-Q16\\convert.exe\" -delay 1x%d %s*.png %s.%s"
# animation::saveGIF({
# best.ml <- optim(c(0, 1), m.ml, data = df) 
# }, movie.name = "animation.gif", interval = 0.5, ani.width = 550, ani.height = 350, convert=imconvertstring)
# 
# ani.options(convert = '/Applications/ImageMagick-7.0.5/bin/convert')
# im.convert("plot_*.png", output="animation.gif")
# 
# plot(df$x, df$y)
# abline(coef=best.ml$par, col = "red")

# we calculate a "classical" lm()
m.lm <- lm(y ~ x, data=df)

print("compare parameters:")
## [1] "compare parameters:"
print("lm parameters: ")
## [1] "lm parameters: "
print(m.lm$coefficients)
## (Intercept)           x 
##    9.856139    1.644558
print("lsq parameters: ")
## [1] "lsq parameters: "
print(best.lsq$par)
## [1] 9.856692 1.644704
print("ml parameters: ")
## [1] "ml parameters: "
print(best.ml$par)
## [1] 9.855540 1.644382

Animation der Anpassung:

least square optimization

least square optimization process

maximum likelihood optimization

maximum likelihood optimization process

Modelltests für das vorgestellte Vorgehen

  • Nullmodel: nur Konstante
  • erweitertes Modell: Konstante und Steigung werden geschätzt

Nicht wirklich sinnvoll, aber zeigt die Universalität des Ansatzes.

Bei einfacher, auch multipler Regression sind die Residuen normalverteilt mit einem Mittelwert von 0.

# generate sample data
set.seed(42)
dd <- rnorm(100, 0, 2)
require(dplyr)
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
df <- dplyr::data_frame(
    y = 10 + dd * 2,
    x = dd + rnorm(100, 0, 1),
    c = rep(1, 100)
)
## Warning: `data_frame()` was deprecated in tibble 1.1.0.
## Please use `tibble()` instead.
# we read a data file created like this to get results comparable
# write.table(ddf, "ml_lm_simple_1.txt", quote=F, sep="\t", row.names=F)
df <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/ml_lm_simple_1.txt")
cor(df$x, df$y)
## [1] 0.8880586
require(psych)
## Loading required package: psych
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
describe(df)
##   vars   n  mean   sd median trimmed  mad   min   max range  skew kurtosis   se
## y    1 100  9.65 4.25   9.32    9.65 3.96 -0.83 19.46 20.29  0.04    -0.22 0.42
## x    2 100 -0.13 2.29  -0.05   -0.08 2.08 -6.57  5.98 12.55 -0.15     0.41 0.23
## c    3 100  1.00 0.00   1.00    1.00 0.00  1.00  1.00  0.00   NaN      NaN 0.00
# first we do a least square linear model
lm.1 <- lm(df$y ~ df$x)
summary(lm.1)
## 
## Call:
## lm(formula = df$y ~ df$x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.0298 -1.2722 -0.0842  1.1775  4.7193 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.8561     0.1966   50.14   <2e-16 ***
## df$x          1.6446     0.0860   19.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.963 on 98 degrees of freedom
## Multiple R-squared:  0.7886, Adjusted R-squared:  0.7865 
## F-statistic: 365.7 on 1 and 98 DF,  p-value: < 2.2e-16
plot(df$x, df$y)
abline(lm.1, col = "red")

#require(tidyverse)
#parms <- tribble( 
#   ~b0, ~m.yHat, ~m.error, 
#   NA, NA, NA
#   )


# we define a function with one parameter to optimize
# we suppose normal distributed error with mean = 0 and put error sd to calculate density for errors
m.b0 <- function(b0) 
    with(df, {
        yHat = b0 + x
        error = y - yHat
        output = -sum(log(dnorm(error, mean = 0, sd=sd(error))))
        #plot(x,y)
        #abline(h=b0, col = "red")
        #print(paste("b0: ", b0, " mean yHat: ", mean(yHat), " error: ", mean(error)))
        #parms <- rbind(parms, c(b0, mean(yHat), mean(error)))
        print(c(b0, mean(yHat), mean(error)))
        return(output)
        }
    )


# we can do MLE using optim(), first parameter is a list of start values for parameters to optimize
# optim(c(0), m.b0, method="BFGS")
#pp <- optim(c(0), m.b0, method="Nelder-Mead")
pp <- optim(c(0), m.b0)
## Warning in optim(c(0), m.b0): one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## [1]  0.000000 -0.125035  9.775547
## [1]  0.10000000 -0.02503496  9.67554681
## [1] 0.20000000 0.07496504 9.57554681
## [1] 0.300000 0.174965 9.475547
## [1] 0.500000 0.374965 9.275547
## [1] 0.700000 0.574965 9.075547
## [1] 1.100000 0.974965 8.675547
## [1] 1.500000 1.374965 8.275547
## [1] 2.300000 2.174965 7.475547
## [1] 3.100000 2.974965 6.675547
## [1] 4.700000 4.574965 5.075547
## [1] 6.300000 6.174965 3.475547
## [1] 9.5000000 9.3749650 0.2755468
## [1] 12.700000 12.574965 -2.924453
## [1] 12.700000 12.574965 -2.924453
## [1] 11.100000 10.974965 -1.324453
## [1] 7.900000 7.774965 1.875547
## [1] 10.3000000 10.1749650 -0.5244532
## [1] 8.700000 8.574965 1.075547
## [1]  9.9000000  9.7749650 -0.1244532
## [1] 10.3000000 10.1749650 -0.5244532
## [1] 9.70000000 9.57496504 0.07554681
## [1] 9.5000000 9.3749650 0.2755468
## [1]  9.80000000  9.67496504 -0.02445319
## [1]  9.9000000  9.7749650 -0.1244532
## [1] 9.75000000 9.62496504 0.02554681
## [1]  9.85000000  9.72496504 -0.07445319
## [1] 9.7750000000 9.6499650362 0.0005468085
## [1] 9.75000000 9.62496504 0.02554681
## [1]  9.78750000  9.66246504 -0.01195319
## [1] 9.76250000 9.63746504 0.01304681
## [1]  9.781250000  9.656215036 -0.005703192
## [1] 9.768750000 9.643715036 0.006796808
## [1]  9.778125000  9.653090036 -0.002578192
## [1] 9.771875000 9.646840036 0.003671808
## [1]  9.776562500  9.651527536 -0.001015692
pp
## $par
## [1] 9.775
## 
## $value
## [1] 230.9785
## 
## $counts
## function gradient 
##       36       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
plot(df$x,df$y)
abline(h=pp$par)

# we can use mle() of package stats4 to do the same and generate a result object
#require(stats4)
require(bbmle)
## Loading required package: bbmle
## 
## Attaching package: 'bbmle'
## The following object is masked from 'package:dplyr':
## 
##     slice
r.b0 = bbmle::mle2(m.b0, start=list(b0=1), control = list(maxit =  10^9, factr = 10^-15))
## [1] 1.000000 0.874965 8.775547
## [1] 1.001000 0.875965 8.774547
## [1] 0.999000 0.873965 8.776547
## [1]  147.2689  147.1438 -137.4933
## [1]  30.25378  30.12874 -20.47823
## [1] 6.850755 6.725720 2.924792
## [1] 6.851755 6.726720 2.923792
## [1] 6.849755 6.724720 2.925792
## [1]  9.775547e+00  9.650512e+00 -4.309304e-12
## [1]  9.776547  9.651512 -0.001000
## [1] 9.774547 9.649512 0.001000
## [1]  9.775547e+00  9.650512e+00 -4.604763e-14
## [1] 9.775547e+00 9.650512e+00 7.100823e-11
## [1] 9.775547e+00 9.650512e+00 7.100823e-11
## [1] 9.775547e+00 9.650512e+00 7.100823e-11
## [1] 10.7531015 10.6280665 -0.9775547
## [1] 8.7979921 8.6729572 0.9775547
## [1] 10.2643241 10.1392892 -0.4887773
## [1] 9.2867695 9.1617345 0.4887773
## [1] 10.0199355  9.8949005 -0.2443887
## [1] 9.5311581 9.4061232 0.2443887
## [1]  9.8977411  9.7727062 -0.1221943
## [1] 9.6533525 9.5283175 0.1221943
## [1] 9.775547e+00 9.650512e+00 7.100823e-11
## [1]  9.7765243631  9.6514893993 -0.0009775546
## [1] 9.7745692537 9.6495342899 0.0009775548
## [1]  9.7760355857  9.6510006219 -0.0004887773
## [1] 9.7750580311 9.6500230673 0.0004887774
## [1]  9.7757911971  9.6507562333 -0.0002443886
## [1] 9.7753024197 9.6502674559 0.0002443887
## [1]  9.7756690027  9.6506340389 -0.0001221943
## [1] 9.7754246141 9.6503896503 0.0001221944
summary(r.b0)
## Maximum likelihood estimation
## 
## Call:
## bbmle::mle2(minuslogl = m.b0, start = list(b0 = 1), control = list(maxit = 10^9, 
##     factr = 10^-15))
## 
## Coefficients:
##    Estimate Std. Error z value     Pr(z)    
## b0  9.77555    0.24494   39.91 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 461.957
r.b0.pred   <- df$x + r.b0@coef[1]
r.b0.err    <- df$y - r.b0.pred
# we estimated the constant only, so our sd of error is almost equal to sd of predictor
sd(r.b0.err)
## [1] 2.449408
sd(df$x)
## [1] 2.293818
AIC(r.b0)
## [1] 463.957
# we define a function with tow parameters constant and slope
m.b1 <- function(b0, b1) with(df, {
        yHat = b0 + b1 * x
        error = y - yHat
        output = -sum(log(dnorm(error, mean = 0, sd=sd(error))))
        return(output)
        }
    )
require(bbmle)
r.b1 = bbmle::mle2(m.b1, start=list(b0=1, b1=1), control = list(maxit =  10^9, factr = 10^-15))
summary(r.b1)
## Maximum likelihood estimation
## 
## Call:
## bbmle::mle2(minuslogl = m.b1, start = list(b0 = 1, b1 = 1), control = list(maxit = 10^9, 
##     factr = 10^-15))
## 
## Coefficients:
##    Estimate Std. Error z value     Pr(z)    
## b0 9.856139   0.195575  50.396 < 2.2e-16 ***
## b1 1.644558   0.085136  19.317 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 416.6461
# and to get the deviance etc
# we want to see the likelihood
exp(logLik(r.b1))
## 'log Lik.' 3.3609e-91 (df=2)
# this is very small and not easy to handle. So we use log(likelihood) usually.
logLik(r.b1)
## 'log Lik.' -208.323 (df=2)
# or even more common -2LL
-2*logLik(r.b1)
## 'log Lik.' 416.6461 (df=2)
# this is, what we call deviance
# in most lm() commands, deviance is a result object attribute
# but this is not true for bbmle result objects, therefore commented out
# r.b1$deviance
# ... and it is in known result objects in R
# sum(residuals(r.b1)^2)

# compare AICs
AIC(r.b0)
## [1] 463.957
AIC(r.b1)
## [1] 420.6461
# we do a model test
anova(r.b0, r.b1)
## Likelihood Ratio Tests
## Model 1: r.b0, [m.b0]: b0
## Model 2: r.b1, [m.b1]: b0+b1
##   Tot Df Deviance  Chisq Df Pr(>Chisq)    
## 1      1   461.96                         
## 2      2   416.65 45.311  1  1.681e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Beispiele, Aufgaben, Übungen einfache Regression

Berechnen Sie die Koeffizienten (Parameter) auf die drei oben vorgestellten Arten und Weisen (lm, least square, maximum likelihood) für die Vorhersage des Körpergewichtes (weight) aufgrund der Körpergröße (height).

df <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_bmi_preprocessed.txt", fileEncoding = "UTF-8")

Likelihood and deviance

In a perfect model the probability of the data with the parameters chosen and their values is 1. This is valid for every single observation. A model, that does not help at all to predict a reaction variable would have a probability of the data given the model of 0. The smaller (worse) the probability the smaller the resulting number. Probabilities vary between 1 and 0. The log of 1 is zero. The log of 0 is \(-inf\). So the smaller (worse) the probability, the more negative is the log of it.

Likelihood is the product of every single likelihood of our observations. log(likelihood) is the sum of the log(likelihood) of our observations.

The bigger the likelihood, the closer the log(likelihood) to 0. The worse the likelihood, the more negative the log(likelihood). -2LL is called the deviance.

Deviance is the sum of the squared residuals in linear models. So the log(likelihood) is the sum of the logs of the squared residuals in a dataset, given the model. So we can consider the deviance to be a measure of how far the data are of beeing perfectly explained by our model.

Therefore the bigger the deviance the less perfect is the model to explain the data. This is valid for the AIC and BIC also, because they depend on the deviance, they modify the deviance.

df <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")
m.mv <- lm(univ_GPA ~ math_SAT + verb_SAT, data=df)
df$pred <- predict(m.mv)
cor(df$univ_GPA, df$pred)
## [1] 0.6857272
deviance(m.mv)
## [1] 11.0184
sum(residuals(m.mv)^2)
## [1] 11.0184
AIC(m.mv)
## [1] 69.26579
BIC(m.mv)
## [1] 79.88163
# comparison 0 model with model chosen
m.0 <- lm(univ_GPA ~ 1, data=df)
# deviance increases
cat(deviance(m.0), '  ', deviance(m.mv))
## 20.79814    11.0184
# AIC decreases
cat(AIC(m.0), '  ', AIC(m.mv))
## 131.9719    69.26579
# BIC decreases
cat(BIC(m.0), '  ', BIC(m.mv))
## 137.2799    79.88163