aspettl — May 10, 2013, 10:37 AM
handmessungen <- read.table("handmessungen.txt", header=T)
# a)
modell <- lm(Groesse ~ Handlaenge + Handbreite, data = handmessungen)
summary(modell)
Call:
lm(formula = Groesse ~ Handlaenge + Handbreite, data = handmessungen)
Residuals:
Min 1Q Median 3Q Max
-15.906 -3.602 -0.683 4.113 19.279
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 69.192 7.681 9.01 1.1e-14 ***
Handlaenge 4.433 0.594 7.46 2.7e-11 ***
Handbreite 1.176 0.440 2.67 0.0087 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.05 on 104 degrees of freedom
Multiple R-squared: 0.647, Adjusted R-squared: 0.641
F-statistic: 95.5 on 2 and 104 DF, p-value: <2e-16
# für später:
X <- model.matrix(modell) # Designmatrix holen/erstellen lassen
Y <- handmessungen$Groesse
n <- length(Y)
beta <- coef(modell) # die geschätzten Regressionskoeffizienten (alternativ modell$coefficients)
m <- length(beta)
S <- summary(modell)$sigma # Standardabweichung der normalverteilten Störgrößen ist 6.049
# alternativ kann man alles auch selbst berechnen:
X <- cbind(rep(1, length(handmessungen$Groesse)), handmessungen$Handlaenge, handmessungen$Handbreite)
Y <- handmessungen$Groesse
n <- length(Y)
beta <- solve(t(X) %*% X) %*% t(X) %*% Y # Schätzer für die Regressionskoeffizienten
m <- length(beta)
S <- sqrt(1/(n-m) * t(Y - X %*% beta) %*% (Y - X %*% beta)) # Schätzer für die Standardabweichung der Störgrößen
# b)
# Nullhypothese H0: Einfluss der Handbreite = 0 (entspricht "nur Handlänge ist relevant")
# Laut summary(modell) ist der Wert der Testgröße 2.675 (Feld "t value") der Handbreite-Zeile.
alpha <- 0.05
qt(1-alpha/2, n-m)
[1] 1.983
# Das (1-alpha/2)-Quantil der t-Verteilung mit n-m Freiheitsgraden ist ca. 1.983.
# Die Nullhypothese wird also abgelehnt, da 2.675 > 1.983.
#
# Man kann auch dies direkt aus summary(modell) ablesen:
# Wenn H0 gilt, dann tritt dieser Wert für die Testgröße mit Wkt. 0.009 (Feld "Pr(>|t|)") auf.
# Da 0.009 < 0.05 ist, wird H0 abgelehnt. Auch bei alpha=0.01 würde noch (d.h. recht knapp) abgelehnt.
#
# alternativ, Berechnung der Testgröße(n):
tstats <- beta / (S * sqrt(diag(solve(t(X) %*% X)))) # hier: Testgrößen für alle drei Komponenten
tstats[3] # Testgröße für Handbreite, also H0
[1] 2.675
# Testgröße (für Handbreite) = 2.675, Vergleich mit Quantil wie oben
# c)
predict(modell, data.frame(Handlaenge = 18.4, Handbreite = 20.3), interval="confidence", level = 0.9)
fit lwr upr
1 174.6 173.7 175.6
# Das 90%-Konfidenzintervall für die mittlere Körpergröße ist [173.7, 175.6].
# alternativ:
x0 <- c(1, 18.4, 20.3) # x-Werte für Vorhersage
z0 <- qt(1-0.1/2, length(handmessungen$Groesse)-3) * S * sqrt(t(x0) %*% solve(t(X) %*% X) %*% x0)
lowerBound = t(beta) %*% x0 - z0
upperBound = t(beta) %*% x0 + z0
print(c(lowerBound, upperBound))
[1] 173.7 175.6