Aflevering 11
Kamille Møller Frost
2024-05-05
Opgave 2 - Reeksamen 2020
Vi betragter i denne opgave laks fanget i Lunefloden i England. For hver laks er køn, logaritmen til længden og logaritmen til vægten registreret i dataen Laks.csv.
Vi antager desuden i opgaven at for en given længde og et givet køn, da er logaritmen til vægten normalfordelt.
Laks <- read.csv("/Users/kamille/Library/CloudStorage/OneDrive-Personligt/Dokumenter/Mat-øk - au/2. år/4. sem - Matematisk statistik/Afleveringer/Afl11-MaSta2020-reeksamen-filer/Laks.csv")
head(Laks)summary(Laks)## Koen LnLaengde LnVaegt
## Length:78 Min. :3.010 Min. :1.100
## Class :character 1st Qu.:3.212 1st Qu.:1.742
## Mode :character Median :3.280 Median :1.950
## Mean :3.271 Mean :1.937
## 3rd Qu.:3.310 3rd Qu.:2.105
## Max. :3.510 Max. :2.770a)
Vi skal nu indføre notation for de målte værdier, og opskrive den statistiske model, hvor hvert køn har sin egen lineære sammenhæng mellem middelværdien af logaritmen til vægten og den forklarende variabel logaritmen til længden, og hvert køn har sin egen varians omkring regressionslinjen.
Herefter skal dataen opdeles efter køn, hvorefter der skal udarbejdes et skøn over variansen omkring regressionslinjen. Til sidst skal der laves et test for, at variansen omkring regressionslinjen er den samme for de to køn.
Vi lader \(Y_i\) betegne logaritmen til vægten for den \(i\)’te laks, vi lader \(x_i\) betegne den logaritmiske længde for den \(i\)’te laks. For at differentiere mellem middelværdierne for de to køn, lader vi \(G_i\) betegne kønnet for den \(i\)’te laks således at \(G_i = 0\) for hunner og \(G_i = 1\) for hanner. Derfor betegner vi \(\mu_{G_i}\) som middelværdien for den logaritmiske vægt for laks af kønnet \(G_i\) for en given længde. Til sidst lader vi \(\sigma_{G_i}^2\) betegne variansen omkring regressionslinjen for laks af kønnet \(G_i\).
Da kan vi opskrive den statistiske model som følgende:
\[Y_i = \mu_{G_i}+\varepsilon_i\] For \(i=1,..., 78\).
Hvor \[\mu_i = \beta_0 + \beta_1 x_i\] Og \(\varepsilon_1, ..., \varepsilon_{78} \; \sim^{i.i.d.} \; N(0, \sigma_{G_i}^2)\) idet hvert køn har sin egen varians omkring linjen.
Vi laver nu anden del af delopgave a), hvor testen skal udføres.
Vi starter med at opdele data efter køn:
hun_laks <- subset(Laks, Koen == "hun")
han_laks <- subset(Laks, Koen == "han")For at lave et skøn over variansen omkring regressionslinjen laver vi lineær regression:
Y_hun <- hun_laks$LnVaegt
Y_han <- han_laks$LnVaegt
x_hun <- hun_laks$LnLaengde
x_han <- han_laks$LnLaengde
M_hun <- lm(Y_hun ~ x_hun, data = hun_laks)
M_han <- lm(Y_han ~ x_han, data = han_laks)
summary(M_hun)##
## Call:
## lm(formula = Y_hun ~ x_hun, data = hun_laks)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.27392 -0.04198 0.00677 0.03965 0.36709
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.2846 0.6635 -12.49 2.98e-14 ***
## x_hun 3.1376 0.2045 15.34 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09782 on 34 degrees of freedom
## Multiple R-squared: 0.8738, Adjusted R-squared: 0.8701
## F-statistic: 235.4 on 1 and 34 DF, p-value: < 2.2e-16summary(M_han)##
## Call:
## lm(formula = Y_han ~ x_han, data = han_laks)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15704 -0.03701 -0.01046 0.02953 0.15439
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.9352 0.5313 -18.70 <2e-16 ***
## x_han 3.6149 0.1612 22.43 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07977 on 40 degrees of freedom
## Multiple R-squared: 0.9263, Adjusted R-squared: 0.9245
## F-statistic: 502.9 on 1 and 40 DF, p-value: < 2.2e-16Herfra ses det, at for hun laksene er Residual standard error 0,09782, som netop er skønnet over variansen omkring regressionslinjen. For han laksene er skønnet over variansen omkring regressionslinjen 0,07977.
Nu mangler jeg blot at lave et test for, at variansn omkring regressionslinjen er den samme for de to køn. Vi vil altså teste følgende nulhypotese:
\[H_0 \; : \qquad \sigma_0^2 = \sigma_1^2\]
For at teste denne nulhypotese benytter vi F-test:
set.seed(42)
res_hun <- M_hun$residuals
res_han <- M_han$residuals
V_hun <- var(res_hun)
V_han <- var(res_han)
df_hun <- length(res_hun) - 1
df_han <- length(res_han) - 1
F_stat <- V_hun / V_han
p_lille <- pf(F_stat, df_hun, df_han)
p_stor <- pf(F_stat, df_hun, df_han, lower.tail = FALSE)
p_value <- 2 * min(p_lille, p_stor)
print(p_value)## [1] 0.2135819Da p-værdien er forholdsvæis høj kan vi ikke forkaste nulhypotesen om at varianserne omkring regressionslinjen for de to køn er ens.
b)
Vi skal nu betragte modellen hvor hvert køn har sin egen lineære sammenhæng mellem middelværdien af logaritmen til vægten og den forklarende variabel logaritmen til længden, og hvor de to køn har samme varians omkring regressionslinjen.
Vi skal med andre ord betragte modellen fra før, hvor vi blot har at
\[\sigma_0^2 = \sigma_1^2\]
Vi starter med at definere den nye statistiske model i R, ved at sætte Y til den logaritmiske vægt og x til den logaritmiske længde. Herefter indføres en ny variabel group, således at vi tager højde for de to grupper.
group <- factor(c(rep("hun", 36), rep("han", 42)))
Y <- Laks$LnVaegt
x <- Laks$LnLaengde
M1 <- lm(Y ~ group + x:group - 1, data = Laks)
res_M1 <- residuals(M1)Vi laver nu et normalt qqplot for residualerne i denne model:
qqnorm(res_M1)
qqline(res_M1)
Det ses tydeligt, at punkterne ligger tæt til linjen med enkelte afvigelser, hvorfor det ligner at residualerne er normalfordelt.
c)
Vi skal nu lave et F-test for hypotesen, at der er samme hældning i regressionslinjenerne for de to køn.
Vi har følgende to modeller: \[Y_{hun} = \beta_{00} + \beta_{01}x_{hun} + \varepsilon_{hun}\] \[Y_{han} = \beta_{10} + \beta_{11} x_{han} + \varepsilon_{han}\] Hvor \(\beta_{01}\) er skæringen for modellen for hunnerne, og \(\beta_{11}\) er skæringen for modellen for hannerne.
Med følgende nulhypotese \[H_0 \; : \qquad \beta_{01} = \beta_{11}\]
Med andre ord vil vi undersøge vha. et F-test om \(M_1\) fra før kan reduceres til \(M_2^\ast\), hvortil vi kalder hypotesen:
\[H_{02}^\ast \; : \qquad \beta_{01} = \beta_{11}\]
Derfor definerer vi \(M_2^\ast\):
M2star <- lm(Y ~ x + group, data = Laks)Således at vi kan modellere reduktionen fra \(M_1\) til \(M_2^\ast\):
anova(M2star, M1)Da F-statistikken ikke tager værdier omkring 1, kan nulhypotesen ikke antages at være sand. Og da F-statistikken heller ikke er større end 1 kan vi ikke forkaste nulhypotesen.
d)
Vi skal nu lave et \(95 \%\) konfidensinterval for forskel mellem han-laks og hun-laks i skæring af regressionslinjen. Derefter skal vi vurdere om det kan antages, at de to skæringer er ens.
Vi har nu følgende statistiske model:
\[M_2^\ast \; : \qquad Y_{G_i} \sim N(\alpha_G + \beta x_{G_i}, \sigma^2)\] Som vi desuden definerede i forrige delopgave.
Vi beregner nu konfidensintervallet ved følgende:
hun_intercept <- coef(M2star)["(Intercept):hun"]
han_intercept <- coef(M2star)["(Intercept):han"]
intercept_diff <- hun_intercept - han_intercept
# Jeg kan ikke finde ud af hvad jeg gør forkert her, men det er som om at min M2star kun spytter et intercept ud...
SE <- sqrt(sum(vcov(M2star)[1,1], vcov(M2star)[2,2]))
alpha <- 0.05
df <- df.residual(M2star)
t_value <- qt(1 - alpha/2, df)
lower_bound <- intercept_diff - t_value * SE
upper_bound <- intercept_diff + t_value * SE
print(lower_bound)## <NA>
## NAprint(upper_bound)## <NA>
## NAe)
Vi skal i den forrige model fra d) lave et \(90 \%\) prædikationsinterval for vægten i pund af en hunlaks med en længde på 27 inch.
Vi starter med at tage logaritmen til 27 inch, så det passer med modellen:
value <- log(27)Og jeg benytter funktionen predict
predicted_weight <- predict(M2star, newdata = data.frame(x = value, group = "hun"), interval = "prediction", level = 0.90)
print(predicted_weight[1])## [1] 1.949153print(predicted_weight[3])## [1] 2.104143Det betyder altså at vi har estimeret vægten af en hunlaks på 27 inch med et \(90 \%\) prædikationsinterval til at være intervallet \([1.9492, 2.1041]\).