Statistik (Fach) / E-TESTs 07 (Lektion)

Berechnungen 07 - 12

Diese Lektion wurde von dollinger3000 erstellt.

Lektion lernen

zurück  |  weiter 1 / 1

• 01 Führe für folgenden Datensatz eine Poisson und eine negative-binomiale Regression durch. Um wieviel ist der AIC im negativ binomialen Modell niedriger? (Bemerkung: Wir vergleichen AIC statt deviance, da im negbin 2 Parameter gefittet werden müssen, im Poisson-Modell nur einer. Dass dann die deviance niedriger ist verwundert zunächst einmal nicht.) x <- 12:1 y <- c(3, 4, 4, 0, 3, 12, 1, 10, 9, 0, 3, 49) > require(MASS)> x1 <- 12:1> y1 <- c(3,4,4,0,3,12,1,10,9,0,3,49)> summary(poisReg<-glm(y1~x1, poisson)) AIC: 133.87> summary(nbReg<-glm.nb(y1~x1))AIC: 76.195> 133.87-76.195 [1] 57.675
• 02 Führe eine Poisson-Regression mit folgenden Daten durch. Berechne dann die Pearson-, deviance und studentised residuals. Für den ersten Datenpunkt, welche Art von Residuen hat den negativsten Wert? x = c(1,2,3,4,5,6,7,8,9,10,11,12) y = c(4, 7, 4, 7, 10, 6, 6, 11, 10, 9, 10, 10) a. deviance residuals b. studentised residuals c. Pearson residuals d. zwei haben den gleichen Wert > x2 = c(1,2,3,4,5,6,7,8,9,10,11,12)> y2 = c(4, 7, 4, 7, 10, 6, 6, 11, 10, 9, 10, 10)> Reg<-glm(y2~x2, poisson)> residuals(Reg, type="pearson") > residuals(Reg, type="deviance")> rstudent(Reg)b. studentised residuals
• 03 Um wieviel nimmt die deviance ab, wenn wir für folgenden Datensatz statt eines linearen ein quadratisches Modell fitten? x <- 1:20 y <- c(3.7, 6.8, 4.2, 4.1, 2.8, 4.1, 3.6, 2.8, 2.6, 2.4, 0.4, -1.1, -1.6, -2, -3.4, -5.3, -4.9, -7.3, -8.7,-11.3) > x <- 1:20> y <- c(3.7,6.8,4.2,4.2,2.8,4.1,3.6,2.8,2.6,2.4,0.4,-1.1,-1.6,-2,-3.4,-5.3,-4.9,-7.3,-8.7,-11.3)> summary(glm(y~x, gaussian)) Residual deviance: 47.666 on 18 degrees of freedom > summary(glm(y~poly(x,2), gaussian)) Residual deviance: 11.99 on 17 degrees of freedomResidual deviance: 47.666- 11.999= 35.676
• 04*** Hier ein Datensatz: bsp <- data.frame("Faktor"=rep(LETTERS[1:4], each=10), "Anzahl"=c(3, 3, 4, 7, 2, 7,7, 5, 5, 1, 4, 4, 7, 5, 8, 6, 7, 13, 5, 8, 4, 1, 2, 0, 1, 1, 0, 1, 4, 1, 3, 4, 3, 2,5, 4, 5, 1, 4, 3)) Rechne ein Poisson-GLM und gib den Schätzwert für den Mittelwert des niedrigsten Levels an (auf 2 Nachkommastellen genau). (Also: GLM fitten, rausbekommen, welcher Level den niedrigsten Wert hat, diesen ausrechnen und auf die response-scale rücktransformieren.) > bsp <- data.frame("Faktor"=rep(LETTERS[1:4], each=10), "Anzahl"=c(3, 3, 4, 7, 2, 7,7, 5, 5, 1, 4, 4, 7, 5, 8, 6, 7, 13, 5, 8, 4, 1, 2, 0, 1, 1, 0, 1, 4, 1, 3, 4, 3, 2, 5, 4, 5, 1, 4, 3))> attach(bsp) > regbsp<-glm(Anzahl~Faktor, poisson)> summary(regbsp)FaktorC -1.0761 0.2990 -3.599 0.000319 ***> bspC<-bsp[Faktor=="C",]> bspC > attach(bspC) > mean(Anzahl)[1] 1.5
• 05 Rechne mit folgenden Daten eine GLM. Gib dann den geschätzen Wert der Steigung (auf der link-scale) auf 2 Nachkommastellen genau an. Jannes <- data.frame("highscore.verbessert"=c("nein", "nein", "ja", "ja", "nein", "nein", "ja", "ja", "ja", "ja"), "stunden.gespielt"=c(0.3, 0.5, 0.5, 0.7, 1, 1, 1.5,2, 2.2, 2.5)) (Hinweis: R macht intern aus "ja" 0 und aus "nein" 1. Eine negative Steigung bedeutet also, dass mit mehr Stunden gespielt auch häufiger der high- score verbessert wurde.) > Jannes <- data.frame("highscore.verbessert"=c("nein", "nein", "ja", "ja", "nein","nein", "ja", "ja", "ja", "ja"), "stunden.gespielt"=c(0.3, 0.5, 0.5, 0.7, 1, 1, 1.5,2, 2.2, 2.5))> attach(Jannes) > plot(highscore.verbessert~stunden.gespielt)> regJannes<-glm(highscore.verbessert~stunden.gespielt, binomial)> summary(regJannes)stunden.gespielt -2.344 1.590 -1.474 0.140Ergebnis: -2.34
• 17*** Die Hirschkühe (engl. "hind") auf Rhum haben auch eine Sozialordnung. Hier die Nachkommenen von 10 Hirschkühen (über einen Zeitraum von 10 Jahren) mit niedrigem Sozialstatus ("omegas"): hinds <- data.frame("name"=c("Lucy", "St. Isabel", "Caro", "Ginger", "Spotty", "Judy", "Pat", "Tiny Joe", "Laura", "Antonia"), "male"=c(2, 1, 4, 4, 5, 2, 4, 3, 3, 2), "female"=c(5, 6, 4, 3, 3, 5, 4, 4, 4, 6)) Welcher Prozentsatz/Bruchteil der Nachkommen ist weiblich? Berechne mit Hilfe eines GLMs auf 2 Nachkommastellen genau. (Hinweis: Rücktransformation auf die response-scale nicht vergessen!) > hinds <- data.frame("name"=c("Lucy", "St. Isabel", "Caro", "Ginger", "Spotty",+ "Judy", "Pat", "Tiny Joe", "Laura", "Antonia"),+ "male"=c(2, 1, 4, 4, 5, 2, 4, 3, 3, 2), "female"=c(5, 6, 4, 3, 3, 5, 4, 4, 4, 6))> attach(hinds)> summary(glm(cbind(female, male)~1, binomial))# auf die Reihenfolge von male/female achten! Estimate Std. Error z value Pr(>|z|)(Intercept) 0.3830 0.2368 1.618 0.106> plogis(0.383) # Wahrscheinlichkeit, dass das Kalb weiblich ist [1] 0.5945965Ergebnis: 0.59 oder 59,45%?
• 07*** Im folgenden Datensatz werden monatliche Todesfälle an Lungenkrankheiten in Großbritannien zwischen 1974 und 1979 angegeben (umgebaut aus dem Datensatz "UKLungDeaths".). Gibt es einen Unterschied zwischen Männern und Frauen? Wieviel Prozent der ursprünglichen Varianz wird durch diesen Faktor erklärt? (R2 in % bitte auf 1 Nachkommastelle genau angeben.) (Nachsatz: Wir korrigieren hier nicht für die absolute Anzahl Männer und Frauen in GB, und wir betreiben auch keine Modelldiagnostik. Es geht hier nur um die technischen Fähigkeiten, eine ANOVA durchzuführen.) lungdeaths <- data.frame("gender"=rep(c("male", "female"), each=72), "deaths"=c(as.vector(mdeaths), as.vector(fdeaths))) > lungdeaths <- data.frame("gender"=rep(c("male", "female"), each=72),"deaths"=c(as.vector(mdeaths), as.vector(fdeaths)))> attach(lungdeaths)> View(lungdeaths) #Anova kann durchgeführt werden, da die Gruppen gleich groß sind (je 72)> summary(lm(deaths~gender)) Multiple R-squared: 0.6685= 66.9 % (pro 100)
• 08 Führe für folgenden Datensatz eine Poisson-Regression durch und bestimme dann mittels des anova- Befehls die Signifikanz des Faktors. Gib den X2-Wert ("Chisq") auf vier Nachkommastellen genau an. bsp <- data.frame("Faktor"=rep(LETTERS[1:4], each=10), "Anzahl"=c(3,3,4,6,2,7,7,5,5,1,4,4,7,5,8,6,7,13,5,8,4,1,2,3,2,4,5,2,4,5,3,4,4,3,5,4,5,3,4,3)) > bsp <- data.frame("Faktor"=rep(LETTERS[1:4], each=10), "Anzahl"=c(3,3,4,6,2,7,7,5,5,1,4,4,7,5,8,6,7,13,5,8,4,1,2,3,2,4,5,2,4,5,3,4,4,3,5,4,5,3,4,3))> attach(bsp)> regbsp<-glm(Anzahl~Faktor, poisson)> summary(regbsp) > anova(regbsp, test="Chisq") Pr(>Chi)0.002036---Ergebnis: 0.0020
• 09 Im Datensatz "OrchardSprays" (laden mittels data(OrchardSprays)) werden 8 Pestizide hinsichtlich ihrer Wirksamkeit verglichen. Führe eine ANOVA durch und daran anschließend einen Tukey honest significant difference post-hoc Test (TukeyHSD). Welchen p-Wert hat der Unterschied zwischen treatment "H" und treatment "E" (auf 2 Nachkommastellen genau)? > data(OrchardSprays)> attach(OrchardSprays)> View(OrchardSprays)                  #um zu sehen, wie die Spalten bzw. die Variablen heißen> aovOS<-aov(decrease~treatment) > TukeyHSD(aovOS)              H-E 27.125 -5.169313 59.41931 0.1621616Ergebnis: 0.16
• 10 Welchen Unterschied kann ich als signifikant nachweisen (mit 99%-iger Sicherheit), wenn ich 144 Plots hälftig auf 2 Kalkdüngung und Kontrolle verteile? Voruntersuchungen haben gezeigt, dass die Standardabweichung meines gemessenen Wurzelwachstums 17.2 g/m2 beträgt. Die Teststärke soll ebenfalls hoch sein, nämlich 0.9. (Das bedeutet, wenn ich einen größeren Unterschied finde, dann ist er signifikant mit p < 0.01, aber wenn der Unterschied kleiner ist, dann bin ich sicher (power > 0.9), dass Kalken keinen Effekt hat.) Bitte Antwort auf 2 Nachkommastellen genau.po > power.t.test(n=72, sd=17.2, sig.level=0.01, power=0.9, delta=NULL) delta = 11.19
• 11 Ich finde bei einer Untersuchung keinen signifikanten Unterschied (0.4 und 1.4, p >0.05) zwischen zwei Behandlungen. Wegen meiner großen Stichprobe (jeweils 1322) ist meine power sehr hoch (0.85). Wie stark ist offensichtlich die VARIANZ (nicht die Standardabweichung!) in meinen Daten? (Achtung! Hier muss das Argument sd=NULL benutzt werden, da sonst die Grundeinstellung sd=1 benutzt wird!! Bitte auf 1 Nachkommastelle genau.) > power.t.test(n=1322, delta=1, sig.level=0.05, power=0.85, sd=NULL) wo-sample t test power calculation n = 1322 delta = 1 sd = 8.577136 sig.level = 0.05 power = 0.85 alternative = two.sided NOTE: n is number in *each* group > 8.577136^2[1] 73.56726 (Varianz=sd hoch2)
• 18*** In einem Kalium(Phosphor)-Düngeexperiment (laden mittels require(MASS); data(npk)) wurde in einem randomised block design experimentiert. Berechne die Signifikanz des K-Effektes auf "yield" mittels ANOVA, zuerst ohne, dann mit dem block-Effekt. Um wieviel wird der P-Wert für K durch den Blockeffekt verbessert? (Auf 3 Nachkommastellen genau, bitte.) > require(MASS)> data(npk)> attach(npk) > summary(aov(yield~K)) P = 0.116 > summary(aov(yield~K + Error(block)))P = 0.0715 > 0.116-0.0715 Ergebnis: 0.045
• 12 Wir erhalten bei einer Analyse folgende Modellausgabe: ....S.27 Dokument Desktop Welche Höhe (hier übrigens in feet angegeben) hat eine Weihrauchkiefer (Pinus taeda, loblolly pine) im Alter von 22 Jahren, wenn sie aus der Samenkollektion 323 stammt? Bitte auf 1 Nachkommastelle genau. = (Intercept) + (Effekt von Alter*Alter) + (Effekt von 323) + (Interaktion von Alter und 323) = (I)+(EA*A)+(E1)+(E2*A) = (-0.660144) + (2.608216*22) + (- 0.218252) + (0.040378*22) = 57.39067 = 57,4
• 13*** Reduziere folgendes Modell so weit wie möglich. Kein Prädiktor ist Teil des Designs, alle sind nur hoffnungsvolle Kandidaten. Wie groß ist die residual deviance (auf 1 Nachkommastelle)? (Hinweis: Den folgenden R-Code im Block nach R kopieren. fm ist das Startmodell, das es zu vereinfachen gilt.) set.seed(14) A <- runif(20, 0, 100) B <- rnorm(20, 15, 5) C <- runif(20, 0, 1) D <- runif(20, 0, 1000) y <- 2+2*B + 10*B*C + rnorm(20, 0, 10) fm <- glm(y ~ (A+B+C+D)^2) anova(fm, test="F") > set.seed(14)> A <- runif(20, 0, 100)> B <- rnorm(20, 15, 5)> C <- runif(20, 0, 1)> D <- runif(20, 0, 1000)> y <- 2+2*B + 10*B*C + rnorm(20, 0, 10)> fm <- glm(y ~ (A+B+C+D)^2)> anova(fm, test="F")…> fm2 <- update(fm, .~. -A:B)> anova(fm2, test="F")…> fm3 <- update(fm2, .~. -A:C)> anova(fm3, test="F")…> fm4 <- update(fm3, .~. -A:D)> anova(fm4, test="F")…> fm5 <- update(fm4, .~. -B:D)> anova(fm5, test="F")…> fm6 <- update(fm5, .~. -C:D)> anova(fm6, test="F")…> fm7 <- update(fm6, .~. -A)> anova(fm7, test="F")…> fm8 <- update(fm7, .~. -D)> anova(fm8, test="F")…> fmversuch <- glm(y ~ (B+C)^2) (#man kann's auch direkt so eingeben, da A und D nicht signifikant sind und daher wegfallen> anova(fmversuch) ... Ergebnis: 1653,0
• 14*** Der Datensatz warpbreaks (laden mittels: data(warpbreaks)) berichtet die Ergebnisse einer Materialprüfung von Wolle. Gemessen wurde, wie häufig ein Webfehler auftrat (breaks) pro yard Garn. Getestet wurden zwei Wollarten (wool, mit Leveln A und B) bei drei angelegten Spannungen (tension, L (low), M (medium), H (high)). Fitte ein GLM für breaks mit quasi-Poisson-Verteilung und Interaktion zwischen wool und tension. Berechne dann die erwartet Anzahl breaks für Wolltyp B (A) bei mittlerer (hoher) Spannung (M). Gib diesen Wert als ganze Zahl gerundet an. > data(warpbreaks)> attach(warpbreaks)> require(lattice)> fmbreaks<-glm(breaks~wool*tension, family=quasipoisson)> anova(fmbreaks) > summary(lm(fmbreaks)) > 44.556-20.000[1] 24.556 = 25 für andere Aufgabe (Wool B, tension M) 44.556-16.333*1-20.556*1-20.000*0+21.111*1+10.556*0 = 28.778 [woolA + woolB + tensionB + wb:tm]
• 15 Es ist wichtig, die Variablen bei einer PCA zu Standardisieren! Diese Aufgabe soll den Effekt verdeutlichen. Führe für folgenden Datensatz zwei Hauptkomponentenanalysen durch, einmal mit Standardisierung, einmal ohne. Wieviel erklärte Varianz erhält die zweite Achse (PC2) durch die Standardisierung HINZU (nicht in Prozent, sondern als Anteil)? Bitte auf 3 Nachkommastellen genau angeben. set.seed(222) A <- runif(20, 0, 100) B <- rnorm(20, 15, 5) C <- runif(20, 0, 1) D <- runif(20, 0, 1000) E <- rnorm(20, 15, 5) F <- rnorm(20, 15, 5) dats <- cbind(A, B, C, D, E, F) set.seed(222)> A <- runif(20, 0, 100)> B <- rnorm(20, 15, 5)> C <- runif(20, 0, 1)> D <- runif(20, 0, 1000)> E <- rnorm(20, 15, 5)> F <- rnorm(20, 15, 5)> dats <- cbind(A, B, C, D, E, F)> View(dats)> pca1<-prcomp(dats)> pca2<-prcomp(dats, scale=T)> summary(pca1) # proportion of variance - pc2 (Textvorgabe) >summary(pca2) # > 0.254-0.02164 [1] 0.23236
• 16*** Reduziere folgendes Modell so weit wie möglich mittels BIC-Kriteriums (Argument k= … in step). Wie groß ist die redisual deviance (auf 1 Nachkommastelle)? (Hinweis: Den folgenden R-Code im Block nach R kopieren. fm ist das Startmodell, das es zu vereinfachen gilt.) set.seed(14) A <- runif(23, 0, 100) B <- rnorm(23, 15, 5) C <- runif(23, 0, 1) D <- runif(23, 0, 1000) y <- 2+2*B + 10*B*C + rnorm(23, 0, 10) fm <- glm(y~ (A+B+C+D)^2) anova(fm, test="F") > set.seed(14)> A <- runif(23, 0, 100)> B <- rnorm(23, 15, 5)> C <- runif(23, 0, 1)> D <- runif(23, 0, 1000)> y <- 2+2*B + 10*B*C + rnorm(23, 0, 10)> fm <- glm(y~ (A+B+C+D)^2)> anova(fm, test=”F”) > fmstep<-step(fm, k=log(23)) Egebnis: 1139.1          #letzte tabelle none / deviance
• 19*** Ein Kollege in Schweden und Du analysieren denselben Datensatz. Auf einem Workshop vergleicht Ihr Eure Ergebnisse. Du hast die Interaktion der zwei Prädiktoren A und B im Modell, er nicht. Um wie viel Prozent unterscheiden sich Eure Modelle für den Punkt (A=10, B=20=? Genauer: Um wie viel Prozent liegt Dein Erwartungswert unter seinem? (Hinweis: Was gesucht ist, ist die Differenz Deines Modells zu seinem relativ zu seinem. Es bedarf folgende Schritte: * Interaktion aus fm entfernen; * mit beiden Modellen auf (10, 20) vorhersagen * Vorhersagewert des Modells ohne Interaktion – Vorhersagewert des Modells mit Interaktion, dann durch Vorhersagewert des Modells ohne Interaktion teilen und mit 100% multiplizieren (oder auch (e-i)/e*100, wobei i der Wert aus dem Interaktionsmodell ist und e aus dem einfachen) * auf 1 Nachkommastellen genau angeben. set.seed (44) A <- round(runif(23, 0, 100), 2) B <- round(rnorm(23, 15, 5), 2) y <- 0.5*A - 0.5*B + 0.25*A*B + rnorm(23, 0,10) fm <- lm(y~A*B) # Modell fitten summary(fm) anova(fm) predict(fm, newdata=data.frame("A"=10, "B"=20)) set.seed (44)A <- round(runif(23, 0, 100), 2)B <- round(rnorm(23, 15, 5), 2)y <- 0.5*A - 0.5*B + 0.25*A*B + rnorm(23, 0,10)fm <- lm(y~A*B) # Modell fittensummary(fm)anova(fm)predict(fm, newdata=data.frame("A"=10, "B"=20)) > fmschwede<-update(fm, .~. -A:B)> summary(fmschwede) > anova(fmschwede) > predict(fmschwede, newdata=data.frame("A"=10, "B"=20))1 73.75381> 73.75381-32.30009[1] 41.45372> (41.45372/73.75381)*100[1] 56.20553 Ergebis: 56.2
• 06 Berechne den Standardfehler des Mittelwerts für den folgenden Datensatz (auf 2 Nachkommastellen genau)! 20.8, 18.3, 23.7, 27.2, 24, 24.4, 18.9, 24.1, 20.2, 26.3, 20.9, 15.8, 18.8 Der Wert muss zwischen 0.8971279 und 0.9971279 liegen > a<-c(20.8, 18.3, 23.7, 27.2, 24, 24.4, 18.9, 24.1, 20.2, 26.3, 20.9, 15.8, 18.8) >install.packages("plotrix") > library(plotrix)> std.error(a)[1] 0.9471279 Der Wert muss zwischen 0.8971279 und 0.9971279 liegen

zurück  |  weiter 1 / 1