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
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
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
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
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.
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
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
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... \]
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
Eine Normalverteilung wird vollständig beschrieben durch Mittelwert und Standardabweichung, also durch die beiden Parameter \(\bar{x}, sd\).
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
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.
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
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
least square optimization
maximum likelihood optimization
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
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")
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