Es gibt einige Datensätze, die in R hinterlegt sind und die man sich zum Üben einfach ins Environment laden kann, wenn man ihren Namen kennt. Hier findet sich eine Übersicht von einigen Datensätzen (wenn man auf den Namen eines Datensatzes klickt, erscheint eine Beschreibung der Variablen).

Die Datensätze lassen sich über folgenden Befehl ins Environment laden und dann ganz normal bearbeiten:

#name <- as.data.frame(name)

Aufgabe 1: Insektensprays

Datensatz

InsectSprays <- as.data.frame(InsectSprays)

Spray = verschiedene Sprays (A-F)

Count = Anzahl von Insekten auf der Haut

Aufgaben:

Folgende Hypothese wurde aufgestellt: Die Sprays unterscheiden sich in ihrer Wirksamkeit.

Eine einfaktorielle Varianzanalyse bietet sich an, da wir eine ungerichtete Hypothese über Gruppenunterschiede testen wollen. (Eine Regression mit Dummykodierung wäre aber auch möglich.)

Für einfaktorielle Varianzanalyse:

H1: mü(k) =/= mü(k’) für mindestens ein paar k,k’; H0: mü(k) = mü(k’) für alle k,k’

(k bezieht sich hier auf die Insektensprays und die Mittelwerte sind die mittlere Anzahl von Insekten auf der Haut nach der Anwendung)

Für Regression:

H1: R2 Modell > R2 Nullmodell; H0: R2 Modell <= R2 Nullmodell

(Hier ist zu beachten, dass die Hypothesen über das Gesamtmodell aufgestellt werden, nicht über einzelne Prädiktoren)

library(Rcmdr)
## Lade nötiges Paket: splines
## Lade nötiges Paket: RcmdrMisc
## Lade nötiges Paket: car
## Lade nötiges Paket: carData
## Lade nötiges Paket: sandwich
## Lade nötiges Paket: effects
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
## Das Commander GUI wird nur bei interaktiven Sessions gestartet
## 
## Attache Paket: 'Rcmdr'
## Das folgende Objekt ist maskiert 'package:base':
## 
##     errorCondition
library(lmtest)
## Lade nötiges Paket: zoo
## 
## Attache Paket: 'zoo'
## Die folgenden Objekte sind maskiert von 'package:base':
## 
##     as.Date, as.Date.numeric
with(InsectSprays, plotMeans(count, spray, error.bars="conf.int", level=0.95, xlab="Art des Sprays", 
  ylab="Anzahl Insekten", connect=TRUE))

Es sieht stark danach aus, dass die Sprays unterschiedlich wirksam sind. Die Mittelwerte von A,B und F sind eindeutig höher als die von C, D und E. Die Konfidenzintervalle der oberen Mittelwerte (A,B,F) beinhalten die unteren Mittelwerte (C,D und E) eindeutig nicht. Weiterhin scheint es, dass C und D sich noch unterscheiden (Mittelwert von D liegt außerhalb des 95% KIs von C und andersrum).

Varianzanalyse

AnovaModel.1 <- aov(count ~ spray, data=InsectSprays)
summary(AnovaModel.1)
##             Df Sum Sq Mean Sq F value Pr(>F)    
## spray        5   2669   533.8    34.7 <2e-16 ***
## Residuals   66   1015    15.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
with(InsectSprays, numSummary(count, groups=spray, statistics=c("mean", "sd")))
##        mean       sd data:n
## A 14.500000 4.719399     12
## B 15.333333 4.271115     12
## C  2.083333 1.975225     12
## D  4.916667 2.503028     12
## E  3.500000 1.732051     12
## F 16.666667 6.213378     12

Regression

summary(lm(count~spray,data=InsectSprays))
## 
## Call:
## lm(formula = count ~ spray, data = InsectSprays)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.333 -1.958 -0.500  1.667  9.333 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  14.5000     1.1322  12.807  < 2e-16 ***
## sprayB        0.8333     1.6011   0.520    0.604    
## sprayC      -12.4167     1.6011  -7.755 7.27e-11 ***
## sprayD       -9.5833     1.6011  -5.985 9.82e-08 ***
## sprayE      -11.0000     1.6011  -6.870 2.75e-09 ***
## sprayF        2.1667     1.6011   1.353    0.181    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.922 on 66 degrees of freedom
## Multiple R-squared:  0.7244, Adjusted R-squared:  0.7036 
## F-statistic:  34.7 on 5 and 66 DF,  p-value: < 2.2e-16

Die Ergebnisse beider Verfahren zeigen, dass sich die Gruppen unterscheiden. Das sieht man in beiden Fällen am signifikanten F-Test (in der Anova für den Faktor Spray, in der Regression für den Test des Gesamtmodells gegen das Nullmodel - beides ist derselbe Test.) Die H1 kann also angenommen und die H0 verworfen werden. (Eins der Verfahren durchzuführen würde natürlich reichen.)

Ein anderes Forschungsteam, das sich mit denselben Daten beschäftigt, hat eine andere Hypothese über die Wirksamkeit. Dieses Team vermutet, dass Spray A weniger wirksam (= höhere Mengen Insekten auf der Haut) sei als alle anderen Sprays.

Hier wäre eine Kontrastanalyse mit Helmert-Kontrasten geeignet, da diese Kontraste es erlauben, eine Gruppe gegen alle anderen zu testen.

H1: D(A gegen B,C,D,E,F) > 0, H0: D(A gegen B,C,D,E,F) <= 0

(Ein Set von Helmert-Kontrasten wird uns natürlich noch weitere Vergleiche zwischen den anderen Gruppen liefern. Über diese wurden hier aber ja keine inhaltlichen Hypothesen aufgestellt.)

Anlegen von Helmert-Kontrasten im Commander und prüfen, ob die Koeffizienten so sind, wie wir sie haben wollen:

contrasts(InsectSprays$spray) <- "contr.helmert"

contrasts(InsectSprays$spray)
##   [,1] [,2] [,3] [,4] [,5]
## A   -1   -1   -1   -1   -1
## B    1   -1   -1   -1   -1
## C    0    2   -1   -1   -1
## D    0    0    3   -1   -1
## E    0    0    0    4   -1
## F    0    0    0    0    5

Die Kontraste sind noch nicht richtig angelegt. Derzeit würde uns der letzte Kontrast (5) die Gruppe F gegen alle anderen Gruppen vergleichen. Wir müssen also die Anordnung der Faktorstufen ändern.

(Commander, Datenmanagement, Neuanordnung der Faktorstufen, umgekehrte Reihenfolge angeben)

InsectSprays$spray <- with(InsectSprays, factor(spray, levels=c('F','E','D','C','B','A')))

Nochmal Kontraste anlegen und prüfen, ob es stimmt:

contrasts(InsectSprays$spray) <- "contr.helmert"

contrasts(InsectSprays$spray)
##   [,1] [,2] [,3] [,4] [,5]
## F   -1   -1   -1   -1   -1
## E    1   -1   -1   -1   -1
## D    0    2   -1   -1   -1
## C    0    0    3   -1   -1
## B    0    0    0    4   -1
## A    0    0    0    0    5

Jetzt stimmt es. Kontrast 5 ist der, auf den sich unsere Hypothese bezieht.

Regression mit kontrastkodierter Variable:

summary(lm(count~spray,data=InsectSprays))
## 
## Call:
## lm(formula = count ~ spray, data = InsectSprays)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.333 -1.958 -0.500  1.667  9.333 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.5000     0.4622  20.554  < 2e-16 ***
## spray1       -6.5833     0.8006  -8.223 1.05e-11 ***
## spray2       -1.7222     0.4622  -3.726 0.000405 ***
## spray3       -1.5694     0.3268  -4.802 9.39e-06 ***
## spray4        1.7083     0.2532   6.748 4.53e-09 ***
## spray5        1.0000     0.2067   4.838 8.22e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.922 on 66 degrees of freedom
## Multiple R-squared:  0.7244, Adjusted R-squared:  0.7036 
## F-statistic:  34.7 on 5 and 66 DF,  p-value: < 2.2e-16

Die H1 kann angenommen werden. Wenn man Spray A gegen alle anderen Gruppen zusammengenommen vergleicht, findet man einen signifikanten Unterschied (Kontrast “spray5” im Output). Nach der Anwendung von Spray A sind also noch mehr Insekten auf der Haut als im Schnitt bei allen anderen Sprays.

Aufgabe 2: Zahnwachstum bei Meerschweinchen

Datensatz

ToothGrowth <- as.data.frame(ToothGrowth)

dose = Dosierung von Vitamic C Supplements, die dem Meerscheinchen gegeben wurden (in milligramm)

supp = Art der Verabreichung (OJ = Orange Juice, VC = Ascorbinsäure)

len = Zahnlänge

Aufgaben:

Es wird vermutet, dass eine höhere Dosierung Vitamin C zu mehr Zahnwachstum (längeren Zähnen) führt. Dabei soll auch für die Art der Verabreichung und mögliche Interaktionen kontrolliert werden.

Am besten passt eine einfache Regression mit Dosierung und Verabreichungsform als Prädiktoren (mit Interaktion). Diese erlaubt eine direkte Testung der gerichteten Hypothese über den Prädiktor Dosierung und kontrolliert gleichzeitig für den Einfluss der Verabreichungsform und eventuelle Interaktionen zwischen Dosierung und Verabreichungsform.

H1: beta(Dosierung) > 0, H0: beta(Dosierung) <= 0

scatterplot(len~dose | supp, regLine=TRUE, smooth=FALSE, boxplots=FALSE, by.groups=TRUE, data=ToothGrowth)

Die Zahnlänge steigt mit höheren Dosen von Vitamin C, wie vermutet. Weiterhin scheint die Gabe über Orangensaft effektiver zu sein. Bei der höchsten Dosis in diesem Design (2 mg) macht die Verabreichungsform offenbar keinen Unterschied mehr. Wir können also erwarten, ein signifikantes und positives beta für Dosis zu finden, ein signifikantes beta für die Art der Verabreichungsform (je nach Reihenfolge der Faktorstufen positiv oder negativ, aber inhaltlich so, dass OJ zu höherer Zahnlänge führt als VC), und vermutlich auch eine Interaktion.

summary(lm(len~dose*supp,data=ToothGrowth))
## 
## Call:
## lm(formula = len ~ dose * supp, data = ToothGrowth)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.2264 -2.8462  0.0504  2.2893  7.9386 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   11.550      1.581   7.304 1.09e-09 ***
## dose           7.811      1.195   6.534 2.03e-08 ***
## suppVC        -8.255      2.236  -3.691 0.000507 ***
## dose:suppVC    3.904      1.691   2.309 0.024631 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.083 on 56 degrees of freedom
## Multiple R-squared:  0.7296, Adjusted R-squared:  0.7151 
## F-statistic: 50.36 on 3 and 56 DF,  p-value: 6.521e-16

Die H1 kann angenommen werden, der Regressionskoeffizienz für Dosis ist positiv und signifikant (p-Wert kann aufgrund der gerichteten Hypothese noch halbiert werden).

Interpretation der Koeffizienten (war nicht Teil der Aufgabe, nur zum Verständnis):

Auf Basis der Ergebnisse der ersten Analyse soll nun zielgerichteter getesten werden, a) ob es einen positiven linearen Trend zwischen Dosis und Zahnlänge gibt, und b) ob die Zahnlänge bei der Gruppe mit der höchsten Dosis höher war als bei den Gruppen darunter, und ob sich diese Gruppen ebenfalls noch unterscheiden (mit 1mg > 0.5mg). Berücksichtigen Sie dabei weiterhin auch die Verabreichungsform und mögliche Interaktionen.

ToothGrowth$dose <- factor(ToothGrowth$dose)
  1. lässt sich über polynomiale Kontraste (Trends) testen.

  2. lässt sich über Helmert-Kontraste testen.

  1. H1: D(linearer Trend für dose) > 0, H0: D(linearer Trend für dose) <= 0

  2. H1,1: D(2mg - 1mg und 0.5 mg) > 0, H0,1: D(2mg - 1mg und 0.5 mg) <= 0 H1,2: D(1mg - 0.5 mg) > 0, H0,2: D(1mg - 0.5 mg) <= 0

  1. Trendkontraste
contrasts(ToothGrowth$dose) <- "contr.poly"

summary(lm(len~dose*supp,data=ToothGrowth))
## 
## Call:
## lm(formula = len ~ dose * supp, data = ToothGrowth)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -8.20  -2.72  -0.27   2.65   8.27 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    20.6633     0.6630  31.166  < 2e-16 ***
## dose.L          9.0722     1.1484   7.900 1.43e-10 ***
## dose.Q         -2.4944     1.1484  -2.172 0.034254 *  
## suppVC         -3.7000     0.9376  -3.946 0.000231 ***
## dose.L:suppVC   3.7689     1.6240   2.321 0.024108 *  
## dose.Q:suppVC   2.7312     1.6240   1.682 0.098394 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.631 on 54 degrees of freedom
## Multiple R-squared:  0.7937, Adjusted R-squared:  0.7746 
## F-statistic: 41.56 on 5 and 54 DF,  p-value: < 2.2e-16

H1 kann angenommen werden, es wurde ein positiver linearer Trend detektiert (signifikantes beta für dose.L). Allerdings kann auch ein negativer quadratischer Trend nicht ausgeschlossen werden (signifikantes beta für dose.Q).

Zusatz: Die Interaktionstermini zeigen uns nun, dass der lineare Trend in der VC-Gruppe stärker war als in OJ (was wir in der ersten Analyse ja auch bereits gesehen hatten)

  1. Helmert-Kontraste

Kontraste definieren, Reihenfolge prüfen

contrasts(ToothGrowth$dose) <- "contr.helmert"

contrasts(ToothGrowth$dose)
##     [,1] [,2]
## 0.5   -1   -1
## 1      1   -1
## 2      0    2

Reihenfolge stimmt.

Lineares Modell mit Kontrasten:

summary(lm(len~dose*supp,data=ToothGrowth))
## 
## Call:
## lm(formula = len ~ dose * supp, data = ToothGrowth)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -8.20  -2.72  -0.27   2.65   8.27 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   20.6633     0.6630  31.166  < 2e-16 ***
## dose1          4.7350     0.8120   5.831 3.18e-07 ***
## dose2          2.6983     0.4688   5.756 4.19e-07 ***
## suppVC        -3.7000     0.9376  -3.946 0.000231 ***
## dose1:suppVC  -0.3400     1.1484  -0.296 0.768308    
## dose2:suppVC   1.8900     0.6630   2.851 0.006166 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.631 on 54 degrees of freedom
## Multiple R-squared:  0.7937, Adjusted R-squared:  0.7746 
## F-statistic: 41.56 on 5 and 54 DF,  p-value: < 2.2e-16

H1,1 und H1,2 können angenommen werden, die betas für beide Kontraste sind signifikant (dose1 und dose2).

Zusatz: Die Interaktionen zeigen uns (wenn signifikant), ob die durch die Kontraste kodierten Gruppenunterschiede für VC größer oder kleiner waren als für OJ. Bei VC war der Abstand zwischen der höchsten Dosis und den beiden niedrigeren Dosen also größer als bei OJ.

Aufgabe 3: Drachen

Der IQ von 480 Drachen wurde bestimmt. Weiterhin wurde ihre Körperlänge gemessen sowie ihr Wohnort (Gebirge und konkreter Ort innerhalb des Gebirges) notiert (Beispiel von https://gkhajduk.github.io/2017-03-09-mixed-models/).

Daten:

dragons <- read.csv("dragons.csv", stringsAsFactors=TRUE)

Hier ist mountainRange bereits ein Faktor. Wenn nicht:

dragons$mountainRange <- factor(dragons$mountainRange)

Lineares Modell

summary(lm(testScore~bodyLength*mountainRange,data=dragons))
## 
## Call:
## lm(formula = testScore ~ bodyLength * mountainRange, data = dragons)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -52.294 -10.148   0.252  10.272  42.526 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                      -44.834322  36.122582  -1.241   0.2152   
## bodyLength                         0.377727   0.200533   1.884   0.0602 . 
## mountainRangeCentral             160.671062  80.791468   1.989   0.0473 * 
## mountainRangeEmmental              6.052942  56.853668   0.106   0.9153   
## mountainRangeJulian              180.105673  58.222030   3.093   0.0021 **
## mountainRangeLigurian            118.255251  57.064330   2.072   0.0388 * 
## mountainRangeMaritime            140.017031  86.032711   1.627   0.1043   
## mountainRangeSarntal             116.391013  53.116571   2.191   0.0289 * 
## mountainRangeSouthern             90.017288  52.224808   1.724   0.0854 . 
## bodyLength:mountainRangeCentral   -0.644224   0.399215  -1.614   0.1073   
## bodyLength:mountainRangeEmmental  -0.005913   0.288749  -0.020   0.9837   
## bodyLength:mountainRangeJulian    -0.680563   0.288654  -2.358   0.0188 * 
## bodyLength:mountainRangeLigurian  -0.530430   0.289988  -1.829   0.0680 . 
## bodyLength:mountainRangeMaritime  -0.487847   0.440039  -1.109   0.2682   
## bodyLength:mountainRangeSarntal   -0.408636   0.278861  -1.465   0.1435   
## bodyLength:mountainRangeSouthern  -0.453093   0.289973  -1.563   0.1188   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.91 on 464 degrees of freedom
## Multiple R-squared:  0.5933, Adjusted R-squared:  0.5801 
## F-statistic: 45.12 on 15 and 464 DF,  p-value: < 2.2e-16

Das Modell mit Gebirge und Körperlänge erklärt ca. 59% der Varianz in den IQ-Werten der Drachen. Körperlänge weist (knapp) keinen signifikanten Zusammenhang mit IQ auf. Drachen in unterschiedlichen Gebirgen scheinen unterschiedlich intelligent zu sein, da die Koeffizienten von manchen der Dummy-Variablen für Gebirge signifikant sind. Körperlänge und Gebirge könnten zudem interagieren, da einer der Interaktionsterme signifikant ist (zumindest in diesem Gebirge ist der Zusammenhang zwischen Körperlänge und Intelligenz also anders im Vergleich zur Referenzgruppe).

Normalverteilung der Residuen

#modell nochmal mit namen anlegen
dragon_model <- lm(testScore~bodyLength*mountainRange,data=dragons)

#regressionsstatistiken hinzufügen
dragons<- within(dragons, {
  fitted.dragon_model <- fitted(dragon_model)
  residuals.dragon_model <- residuals(dragon_model)
  rstudent.dragon_model <- rstudent(dragon_model)
  hatvalues.dragon_model <- hatvalues(dragon_model)
  cooks.distance.dragon_model <- cooks.distance(dragon_model) 
})

#histogramm der residuen
with(dragons, Hist(residuals.dragon_model, scale="frequency", breaks="Sturges", col="darkgray"))

#qqplot
with(dragons, qqPlot(residuals.dragon_model, dist="norm", id=list(method="y", n=2, labels=rownames(dragons))))

## [1] 459 148
#shapiro wilk test auf normalverteilung
normalityTest(~residuals.dragon_model, test="shapiro.test", data=dragons)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals.dragon_model
## W = 0.99876, p-value = 0.9854

Residuen sind normalverteilt.

Linearität

statistische Prüfung

resettest(testScore ~ bodyLength * mountainRange, power=2:3, type="regressor", data=dragons)
## 
##  RESET test
## 
## data:  testScore ~ bodyLength * mountainRange
## RESET = 0.20365, df1 = 2, df2 = 462, p-value = 0.8158

nicht signifikant, Linearität kann also angenommen werden

Homoskedastizität

Grafisch:

scatterplot(residuals.dragon_model~fitted.dragon_model, regLine=TRUE, smooth=list(span=0.5, spread=FALSE), boxplots=FALSE, 
  data=dragons)

Statistisch:

bptest(testScore ~ bodyLength * mountainRange, varformula = ~ fitted.values(dragon_model), studentize=FALSE, data=dragons)
## 
##  Breusch-Pagan test
## 
## data:  testScore ~ bodyLength * mountainRange
## BP = 3.6303, df = 1, p-value = 0.05674

Knapp nicht signifikant, Homoskedastizität kann also noch angenommen werden

Verzerrende Ausreißer:

with(dragons, Hist(cooks.distance.dragon_model, scale="frequency", breaks="Sturges", col="darkgray"))

Es scheinen keine verzerrenden Ausreißer vorzuliegen (alle CD-Werte sogar deutlich unter 1)

Führen Sie die Analyse erneut durch, aber diesmal als Multilevel-Modell. Prüfen Sie, ob ein random intercept für Gebirge notwendig ist und/oder ein ramdom slope. Bestimmen Sie das Modell, das die Daten am besten beschreibt. Bestimmen Sie die R2-Annäherung für dieses Modell.

Erstellen Sie vorher auch eine geeignete Grafik und erläutern Siem welche Effekte Sie erwarten.

Benötigte Pakete laden

library(nlme) #für Berechnen des Modells
library(lmtest) #für likelihood ratio tests
library(rcompanion) #für Effektgrößen

Grafik

scatterplot(testScore~bodyLength | mountainRange, regLine=TRUE, smooth=FALSE, boxplots=FALSE, by.groups=TRUE, data=dragons)
## Warning in scatterplot.default(X[, 2], X[, 1], groups = X[, 3], xlab = xlab, : number of groups exceeds number of available colors
##   colors are recycled

Die Intelligenz der Drachen scheint sich nach Gebirge zu unterscheiden, ein random intercept für Gebirge ist vermutlich eine gute Wahl. Die Steigung des IQs mit Körperlänge scheint in manchen Gebirgen positiv (zb “Emmental”) zu sein, in anderen negativ (zb “Julian”). Möglicherweise ist also auch ein random slope für Gebirge gut.

Nullmodell:

Null <- gls(testScore ~ 1, data = dragons, method="ML")

summary(Null)
## Generalized least squares fit by maximum likelihood
##   Model: testScore ~ 1 
##   Data: dragons 
##        AIC      BIC    logLik
##   4375.622 4383.969 -2185.811
## 
## Coefficients:
##                Value Std.Error t-value p-value
## (Intercept) 50.38603  1.050204 47.9774       0
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.19214319 -0.77364791  0.02672957  0.74522002  2.15855302 
## 
## Residual standard error: 22.98483 
## Degrees of freedom: 480 total; 479 residual

Modell mit random intercept für Gebirge, Vergleich mit Nullmodell

RI <- lme(testScore ~ 1, data= dragons, random = ~1|mountainRange, method="ML")
summary(RI)
## Linear mixed-effects model fit by maximum likelihood
##   Data: dragons 
##        AIC      BIC    logLik
##   3999.685 4012.206 -1996.843
## 
## Random effects:
##  Formula: ~1 | mountainRange
##         (Intercept) Residual
## StdDev:    17.46335 14.94435
## 
## Fixed effects:  testScore ~ 1 
##                Value Std.Error  DF  t-value p-value
## (Intercept) 50.38603  6.218273 472 8.102898       0
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.48568314 -0.65893119  0.02154124  0.66645926  2.97295416 
## 
## Number of Observations: 480
## Number of Groups: 8
lrtest(Null,RI)
## Warning in modelUpdate(objects[[i - 1]], objects[[i]]): original model was of
## class "gls", updated model is of class "lme"
## Likelihood ratio test
## 
## Model 1: testScore ~ 1
## Model 2: testScore ~ 1
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   2 -2185.8                         
## 2   3 -1996.8  1 377.94  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modell mit RI ist besser

Modell mit RI und fixed effect Körperlänge, Test gegen Modell mit nur RI

RI_pred <- lme(testScore ~ bodyLength, data= dragons, random = ~1|mountainRange, method="ML")
summary(RI_pred)
## Linear mixed-effects model fit by maximum likelihood
##   Data: dragons 
##        AIC      BIC    logLik
##   4001.478 4018.173 -1996.739
## 
## Random effects:
##  Formula: ~1 | mountainRange
##         (Intercept) Residual
## StdDev:    17.16893 14.94532
## 
## Fixed effects:  testScore ~ bodyLength 
##                Value Std.Error  DF   t-value p-value
## (Intercept) 43.11736 16.961104 471 2.5421317  0.0113
## bodyLength   0.03611  0.078573 471 0.4595173  0.6461
##  Correlation: 
##            (Intr)
## bodyLength -0.933
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.483284485 -0.651994065  0.005146779  0.669704260  2.958968835 
## 
## Number of Observations: 480
## Number of Groups: 8
lrtest(RI,RI_pred)
## Likelihood ratio test
## 
## Model 1: testScore ~ 1
## Model 2: testScore ~ bodyLength
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   3 -1996.8                     
## 2   4 -1996.7  1 0.2075     0.6488

Nicht signifikant besser, fixed effect für Körperlänge fliegt raus.

Modell mit zusätzlich random slope für Gebirge, Test gegen Modell mit nur RI

RI_RS <- lme(testScore ~ 1, data = dragons, random = ~bodyLength|mountainRange, method="ML")

summary(RI_RS)
## Linear mixed-effects model fit by maximum likelihood
##   Data: dragons 
##        AIC      BIC    logLik
##   4003.685 4024.554 -1996.843
## 
## Random effects:
##  Formula: ~bodyLength | mountainRange
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev       Corr  
## (Intercept) 1.746335e+01 (Intr)
## bodyLength  2.025726e-05 0     
## Residual    1.494435e+01       
## 
## Fixed effects:  testScore ~ 1 
##                Value Std.Error  DF  t-value p-value
## (Intercept) 50.38603  6.218274 472 8.102897       0
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.48568315 -0.65893119  0.02154126  0.66645927  2.97295417 
## 
## Number of Observations: 480
## Number of Groups: 8
lrtest(RI_RS,RI)
## Likelihood ratio test
## 
## Model 1: testScore ~ 1
## Model 2: testScore ~ 1
##   #Df  LogLik Df Chisq Pr(>Chisq)
## 1   5 -1996.8                    
## 2   3 -1996.8 -2     0          1

Nicht signifikant besser. Das beste Modell ist also ein Modell mit nur einem random intercept für Gebirge. Die partiellen Unterschiede im Zusammenhang von Körperlänge mit IQ je nach Gebirgeart, die wir in der ersten Analyse beobachtet hatten (die teilweise signifikanten Interaktionsterme) und die wir nach Inspektion des Graphen hätten vermuten können, waren offenbar nicht stark genug, um hier in random slope für Gebirge signifikant werden zu lassen.

R-Quadrat für das finale Modell

nagelkerke(RI,Null)
## $Models
##                                                                     
## Model: "lme.formula, testScore ~ 1, dragons, ~1 | mountainRange, ML"
## Null:  "gls, testScore ~ 1, dragons, ML"                            
## 
## $Pseudo.R.squared.for.model.vs.null
##                              Pseudo.R.squared
## McFadden                            0.0864523
## Cox and Snell (ML)                  0.5449590
## Nagelkerke (Cragg and Uhler)        0.5450200
## 
## $Likelihood.ratio.test
##  Df.diff LogLik.diff  Chisq    p.value
##       -1     -188.97 377.94 3.5007e-84
## 
## $Number.of.observations
##           
## Model: 480
## Null:  480
## 
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
## 
## $Warnings
## [1] "None"

Das finale Modell erklärt ca. 54% der Varianz in den IQ-Werten (Cox and Snell).