10  Wprowadzenie do Analizy Korelacji i Regresji

10.1 Wprowadzenie

Różnica między korelacją (correlation) a przyczynowością/kausalnością (causation) to jedno z podstawowych wyzwań w analizie statystycznej. Korelacja mierzy statystyczny związek między zmiennymi, natomiast przyczynowość oznacza bezpośredni wpływ jednej zmiennej na drugą.

Zależności statystyczne stanowią fundament podejmowania decyzji opartych na danych w wielu dyscyplinach — od ekonomii i zdrowia publicznego po psychologię i nauki o środowisku. Zrozumienie, kiedy związek wskazuje jedynie na asocjację (association), a kiedy na prawdziwą kausalność (genuine causality), jest kluczowe dla poprawnych wniosków i skutecznych rekomendacji politycznych.

10.2 Kowariancja (Covariance)

Kowariancja (covariance) mierzy, w jaki sposób dwie zmienne współzmieniają się, wskazując zarówno kierunek, jak i siłę ich liniowego związku.

Wzór (z próby):

\operatorname{cov}(X,Y) = \frac{\sum_{i=1}^{n} (x_i - \bar x)(y_i - \bar y)}{n - 1}.

Gdzie:

• x_i i y_i to poszczególne obserwacje,
• \bar{x} i \bar{y} to średnie odpowiednio zmiennych X i Y,
• n to liczba obserwacji,
• dzielimy przez (n-1), ponieważ liczymy kowariancję z próby (tzw. poprawka Bessela; Bessel’s correction).

Obliczenia ręczne krok po kroku (Step-by-Step Manual Calculation Process)

Przykład 1: Godziny nauki a wyniki testu

Dane:

• Godziny nauki (X): 2, 4, 6, 8, 10
• Wyniki testu (Y): 65, 70, 80, 85, 95

Krok	Opis	Obliczenia
1	Oblicz średnie	\bar{x}=\frac{2+4+6+8+10}{5}=6 godz.
		\bar{y}=\frac{65+70+80+85+95}{5}=79 pkt
2	Odchylenia od średnich	(x_i-\bar{x}): -4, -2, 0, 2, 4
		(y_i-\bar{y}): -14, -9, 1, 6, 16
3	Iloczyny odchyleń	(x_i-\bar{x})(y_i-\bar{y}): 56, 18, 0, 12, 64
4	Suma iloczynów	\sum = 56+18+0+12+64=150
5	Podziel przez (n-1)	\operatorname{cov}(X,Y)=\frac{150}{5-1}=\frac{150}{4}=37.5

Weryfikacja w R (R Verification):

# Definicja danych
study_hours <- c(2, 4, 6, 8, 10)
test_scores <- c(65, 70, 80, 85, 95)

# Kowariancja
cov_manual <- cov(study_hours, test_scores)
print(paste("Covariance:", round(cov_manual, 2)))

[1] "Covariance: 37.5"

# Weryfikacja krok po kroku
x_mean <- mean(study_hours)
y_mean <- mean(test_scores)
x_dev <- study_hours - x_mean
y_dev <- test_scores - y_mean
products <- x_dev * y_dev
cov_calculated <- sum(products) / (length(study_hours) - 1)

# Tabela obliczeń
data.frame(
  X = study_hours,
  Y = test_scores,
  X_deviation = x_dev,
  Y_deviation = y_dev,
  Product = products
)

   X  Y X_deviation Y_deviation Product
1  2 65          -4         -14      56
2  4 70          -2          -9      18
3  6 80           0           1       0
4  8 85           2           6      12
5 10 95           4          16      64

Interpretacja: Dodatnia kowariancja (37.5) wskazuje, że wraz ze wzrostem liczby godzin nauki rosną także wyniki testu — zmienne mają tendencję do wspólnego wzrostu.

Zadanie ćwiczeniowe z rozwiązaniem (Practice Problem with Solution)

Policz ręcznie kowariancję dla:

• Temperatura (°F): 32, 50, 68, 86, 95
• Sprzedaż lodów ($): 100, 200, 400, 600, 800

Rozwiązanie:

Krok	Obliczenie
1. Średnie	\bar{x}=\frac{32+50+68+86+95}{5}=66.2^{\circ}\mathrm{F}
	\bar{y}=\frac{100+200+400+600+800}{5}=420
2. Odchylenia	X: -34.2, -16.2, 1.8, 19.8, 28.8
	Y: -320, -220, -20, 180, 380
3. Iloczyny	10944, 3564, -36, 3564, 10944
4. Suma	28980
5. Kowariancja	\frac{28980}{4}=7245

# Weryfikacja zadania
temp <- c(32, 50, 68, 86, 95)
sales <- c(100, 200, 400, 600, 800)
print(paste("Covariance for practice problem:", round(cov(temp, sales), 2)))

[1] "Covariance for practice problem: 7245"

10.3 Współczynnik korelacji (Correlation Coefficient)

Współczynnik korelacji (correlation coefficient) standaryzuje kowariancję, usuwając zależność od skali i przyjmując wartości od -1 do +1.

Wskazówki interpretacyjne (Interpretation Guidelines)

Wartość korelacji	Siła	Interpretacja	Przykład
±0.90 do ±1.00	Bardzo silna	Niemal doskonały związek	Wzrost rodziców i dzieci
±0.70 do ±0.89	Silna	Zmienne silnie powiązane	Czas nauki i oceny
±0.50 do ±0.69	Umiarkowana	Umiarkowany związek	Ćwiczenia a spadek masy
±0.30 do ±0.49	Słaba	Słaby związek	Rozmiar buta a umiejętność czytania
±0.00 do ±0.29	Bardzo słaba/brak	Znikomy lub brak związku	Miesiąc urodzenia a inteligencja

Wizualizacja typów zależności korelacyjnych (Types of Correlations Visualization)

# Generowanie przykładowych danych dla różnych wzorców korelacji
n <- 100

# Dodatnia zależność liniowa
set.seed(123)
years_education <- rnorm(n, 14, 3)
annual_income <- 15000 + 3500 * years_education + rnorm(n, 0, 10000)
pos_cor <- round(cor(years_education, annual_income), 3)

# Ujemna zależność liniowa
screen_time <- runif(n, 1, 8)
sleep_hours <- 9 - 0.5 * screen_time + rnorm(n, 0, 0.5)
neg_cor <- round(cor(screen_time, sleep_hours), 3)

# Brak korelacji
x_random <- rnorm(n, 50, 10)
y_random <- rnorm(n, 50, 10)
no_cor <- round(cor(x_random, y_random), 3)

# Nieliniowa zależność (kwadratowa)
hours_studied <- runif(n, 0, 12)
test_score <- -2 * (hours_studied - 6)^2 + 90 + rnorm(n, 0, 5)
nonlin_cor <- round(cor(hours_studied, test_score), 3)

# Ramki danych z wartościami korelacji
positive_data <- data.frame(
  x = years_education, 
  y = annual_income,
  label = paste0("Dodatnia liniowa (r = ", pos_cor, ")")
)

negative_data <- data.frame(
  x = screen_time, 
  y = sleep_hours,
  label = paste0("Ujemna liniowa (r = ", neg_cor, ")")
)

no_corr_data <- data.frame(
  x = x_random, 
  y = y_random,
  label = paste0("Brak korelacji (r = ", no_cor, ")")
)

nonlinear_data <- data.frame(
  x = hours_studied, 
  y = test_score,
  label = paste0("Nieliniowa (Pearson r = ", nonlin_cor, ")")
)

# Połączenie danych
all_data <- rbind(positive_data, negative_data, no_corr_data, nonlinear_data)

# Wykres fasetowy
ggplot(all_data, aes(x = x, y = y)) +
  geom_point(alpha = 0.6, color = "darkblue", size = 2) +
  geom_smooth(method = "lm", color = "red", se = TRUE, alpha = 0.3) +
  facet_wrap(~ label, scales = "free", ncol = 2) +
  labs(
    title = "Różne typy korelacji",
    subtitle = "Linia regresji liniowej (na czerwono) z pasmem ufności",
    x = "", 
    y = ""
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(size = 10, face = "bold"),
    plot.title = element_text(size = 14, face = "bold"),
    plot.subtitle = element_text(size = 10, color = "gray50")
  )

`geom_smooth()` using formula = 'y ~ x'

10.4 Korelacja Pearsona (Pearson Correlation)

Wzór:

r = \frac{\operatorname{cov}(X,Y)}{s_X\, s_Y} = \frac{\sum_{i=1}^n(x_i - \bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i=1}^n(x_i - \bar{x})^2}\,\sqrt{\sum_{i=1}^n(y_i - \bar{y})^2}}.

Pełny przykład obliczeń ręcznych (Complete Manual Calculation Example)

Na danych o godzinach nauki:

• Godziny nauki (X): 2, 4, 6, 8, 10
• Wyniki testu (Y): 65, 70, 80, 85, 95

Szczegółowe kroki:

Krok	Opis	Obliczenia
1	Kowariancja	Z powyżej: \operatorname{cov}(X,Y) = 37.5
2	Kwadraty odchyleń
	Dla X	(x_i-\bar{x})^2: 16, 4, 0, 4, 16
		Suma = 40
	Dla Y	(y_i-\bar{y})^2: 196, 81, 1, 36, 256
		Suma = 570
3	Odchylenia standardowe (standard deviations)
	s_X	s_X=\sqrt{\frac{40}{4}}=\sqrt{10}=3.162
	s_Y	s_Y=\sqrt{\frac{570}{4}}=\sqrt{142.5}=11.937
4	Korelacja	r=\frac{37.5}{3.162 \times 11.937}
		r=\frac{37.5}{37.73}\approx 0.994

# Weryfikacja obliczeń
study_hours <- c(2, 4, 6, 8, 10)
test_scores <- c(65, 70, 80, 85, 95)

# Współczynnik korelacji Pearsona
r_value <- cor(study_hours, test_scores, method = "pearson")
print(paste("Pearson correlation coefficient:", round(r_value, 3)))

[1] "Pearson correlation coefficient: 0.993"

# Obliczenia szczegółowe
x_mean <- mean(study_hours)
y_mean <- mean(test_scores)
x_dev <- study_hours - x_mean
y_dev <- test_scores - y_mean

# Tabela obliczeń
calc_table <- data.frame(
  X = study_hours,
  Y = test_scores,
  X_dev = x_dev,
  Y_dev = y_dev,
  X_dev_sq = x_dev^2,
  Y_dev_sq = y_dev^2,
  Product = x_dev * y_dev
)
print(calc_table)

   X  Y X_dev Y_dev X_dev_sq Y_dev_sq Product
1  2 65    -4   -14       16      196      56
2  4 70    -2    -9        4       81      18
3  6 80     0     1        0        1       0
4  8 85     2     6        4       36      12
5 10 95     4    16       16      256      64

# Statystyki podsumowujące
cat("\nSummary:")

Summary:

cat("\nSum of (X-mean)^2:", sum(x_dev^2))

Sum of (X-mean)^2: 40

cat("\nSum of (Y-mean)^2:", sum(y_dev^2))

Sum of (Y-mean)^2: 570

cat("\nSum of products:", sum(x_dev * y_dev))

Sum of products: 150

cat("\nStandard deviation of X:", sd(study_hours))

Standard deviation of X: 3.162278

cat("\nStandard deviation of Y:", sd(test_scores))

Standard deviation of Y: 11.93734

# Przedział ufności i p-value
cor_test <- cor.test(study_hours, test_scores, method = "pearson")
print(cor_test)

    Pearson's product-moment correlation

data:  study_hours and test_scores
t = 15, df = 3, p-value = 0.0006431
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8994446 0.9995859
sample estimates:
      cor 
0.9933993 

Interpretacja: r \approx 0.994 wskazuje na niemal doskonały dodatni liniowy związek między godzinami nauki a wynikiem testu. Wartość p < 0.05 sugeruje statystyczną istotność tej zależności.

10.5 Korelacja rang Spearmana (Spearman Rank Correlation)

Korelacja Spearmana mierzy monotoniczne zależności, używając rang zamiast surowych wartości.

Wzór:

\rho = 1 - \frac{6 \sum_{i=1}^n d_i^2}{n(n^2 - 1)},

gdzie d_i to różnica rang dla obserwacji i.

Pełny przykład obliczeń ręcznych (Complete Manual Example)

Dane: Wyniki z matematyki i angielskiego

Uczeń	Matematyka	Angielski
A	78	82
B	85	90
C	88	88
D	92	94
E	95	96

Rangowanie i obliczenia:

Uczeń	Wynik z mat.	Ranga mat.	Wynik z ang.	Ranga ang.	d = \text{ranga mat.} - \text{ranga ang.}	d^2
A	78	1	82	1	0	0
B	85	2	90	3	-1	1
C	88	3	88	2	1	1
D	92	4	94	4	0	0
E	95	5	96	5	0	0
					Suma:	2

Obliczenie:

\rho = 1 - \frac{6 \times 2}{5(25-1)} = 1 - \frac{12}{120} = 0.9.

# Dane
math_scores <- c(78, 85, 88, 92, 95)
english_scores <- c(82, 90, 88, 94, 96)

# Rangi
rank_table <- data.frame(
  Student = LETTERS[1:5],
  Math_Score = math_scores,
  Math_Rank = rank(math_scores),
  English_Score = english_scores,
  English_Rank = rank(english_scores)
)

rank_table$d <- rank_table$Math_Rank - rank_table$English_Rank
rank_table$d_squared <- rank_table$d^2

print(rank_table)

  Student Math_Score Math_Rank English_Score English_Rank  d d_squared
1       A         78         1            82            1  0         0
2       B         85         2            90            3 -1         1
3       C         88         3            88            2  1         1
4       D         92         4            94            4  0         0
5       E         95         5            96            5  0         0

cat("\nSum of d^2:", sum(rank_table$d_squared))

Sum of d^2: 2

# Korelacja Spearmana
rho <- cor(math_scores, english_scores, method = "spearman")
cat("\nSpearman correlation (R):", round(rho, 3))

Spearman correlation (R): 0.9

# Obliczenie ręczne
n <- length(math_scores)
rho_manual <- 1 - (6 * sum(rank_table$d_squared)) / (n * (n^2 - 1))
cat("\nSpearman correlation (manual):", round(rho_manual, 3))

Spearman correlation (manual): 0.9

10.6 Tabele krzyżowe (Cross-tabulation) i dane kategoryczne

Tabela krzyżowa (cross-tabulation, contingency table) pokazuje zależności między zmiennymi kategorycznymi.

# Bardziej realistyczne dane przykładowe
set.seed(123)
n_total <- 120

# Poziom edukacji a zatrudnienie
education <- factor(
  c(rep("High School", 40), 
    rep("College", 50), 
    rep("Graduate", 30)),
  levels = c("High School", "College", "Graduate")
)

# Status zatrudnienia z prawdopodobieństwami zależnymi od edukacji
employment <- factor(
  c(
    # High School: 75% employed
    sample(c(rep("Employed", 30), rep("Unemployed", 10)), 40),
    # College: 90% employed  
    sample(c(rep("Employed", 45), rep("Unemployed", 5)), 50),
    # Graduate: 95% employed
    sample(c(rep("Employed", 28), rep("Unemployed", 2)), 30)
  ),
  levels = c("Employed", "Unemployed")
)

# Tabela kontyngencji
contingency_table <- table(education, employment)
print("Contingency Table:")

[1] "Contingency Table:"

print(contingency_table)

             employment
education     Employed Unemployed
  High School       30         10
  College           45          5
  Graduate          28          2

# Procenty w wierszach
print("\nRow Percentages (% within each education level):")

[1] "\nRow Percentages (% within each education level):"

print(round(prop.table(contingency_table, 1) * 100, 1))

             employment
education     Employed Unemployed
  High School     75.0       25.0
  College         90.0       10.0
  Graduate        93.3        6.7

# Test niezależności chi-kwadrat (Chi-square test)
chi_test <- chisq.test(contingency_table)

Warning in chisq.test(contingency_table): Chi-squared approximation may be
incorrect

print(chi_test)

    Pearson's Chi-squared test

data:  contingency_table
X-squared = 5.9623, df = 2, p-value = 0.05073

10.7 Ćwiczenia praktyczne z rozwiązaniami (Practical Exercises with Solutions)

Ćwiczenie 1: Ręczne obliczenie korelacji Pearsona (Calculate Pearson Correlation Manually)

Dane:

• Wzrost (cale): 66, 68, 70, 72, 74
• Waga (funty): 140, 155, 170, 185, 200

Rozwiązanie:

# Dane
height <- c(66, 68, 70, 72, 74)
weight <- c(140, 155, 170, 185, 200)

# Krok 1: Średnie
h_mean <- mean(height)  # 70
w_mean <- mean(weight)  # 170

# Krok 2: Odchylenia i iloczyny
calculation_df <- data.frame(
  Height = height,
  Weight = weight,
  H_dev = height - h_mean,
  W_dev = weight - w_mean,
  H_dev_sq = (height - h_mean)^2,
  W_dev_sq = (weight - w_mean)^2,
  Product = (height - h_mean) * (weight - w_mean)
)

print(calculation_df)

  Height Weight H_dev W_dev H_dev_sq W_dev_sq Product
1     66    140    -4   -30       16      900     120
2     68    155    -2   -15        4      225      30
3     70    170     0     0        0        0       0
4     72    185     2    15        4      225      30
5     74    200     4    30       16      900     120

# Krok 3: Korelacja
cov_hw <- sum(calculation_df$Product) / (length(height) - 1)
sd_h <- sqrt(sum(calculation_df$H_dev_sq) / (length(height) - 1))
sd_w <- sqrt(sum(calculation_df$W_dev_sq) / (length(weight) - 1))
r_manual <- cov_hw / (sd_h * sd_w)

cat("\nManual calculation:")

Manual calculation:

cat("\nCovariance:", round(cov_hw, 2))

Covariance: 75

cat("\nSD Height:", round(sd_h, 2))

SD Height: 3.16

cat("\nSD Weight:", round(sd_w, 2))

SD Weight: 23.72

cat("\nCorrelation:", round(r_manual, 3))

Correlation: 1

# Weryfikacja w R
cat("\n\nR calculation:", round(cor(height, weight), 3))

R calculation: 1

Ćwiczenie 2: Ręczne obliczenie korelacji Spearmana (Calculate Spearman Correlation Manually)

Dane:

• Rangi uczniów z matematyki: 1, 3, 2, 5, 4
• Rangi uczniów z nauk ścisłych (science): 2, 4, 1, 5, 3

Rozwiązanie:

# Rangi (już zrankowane)
math_rank <- c(1, 3, 2, 5, 4)
science_rank <- c(2, 4, 1, 5, 3)

# Różnice
d <- math_rank - science_rank
d_squared <- d^2

# Tabela
rank_df <- data.frame(
  Student = 1:5,
  Math_Rank = math_rank,
  Science_Rank = science_rank,
  d = d,
  d_squared = d_squared
)

print(rank_df)

  Student Math_Rank Science_Rank  d d_squared
1       1         1            2 -1         1
2       2         3            4 -1         1
3       3         2            1  1         1
4       4         5            5  0         0
5       5         4            3  1         1

# Korelacja Spearmana
n <- length(math_rank)
sum_d_sq <- sum(d_squared)
rho_manual <- 1 - (6 * sum_d_sq) / (n * (n^2 - 1))

cat("\nSum of d^2:", sum_d_sq)

Sum of d^2: 4

cat("\nSpearman correlation (manual):", round(rho_manual, 3))

Spearman correlation (manual): 0.8

cat("\nSpearman correlation (R):", round(cor(math_rank, science_rank, method = "spearman"), 3))

Spearman correlation (R): 0.8

Ćwiczenie 3: Interpretacja wyników (Interpretation Practice)

Zinterpretuj następujące wartości korelacji:

1. r = 0.85 między godzinami treningu a wynikiem sprawdzianu sprawności

   Odpowiedź: Silny dodatni związek. Wraz ze wzrostem liczby godzin treningu wyniki istotnie rosną.

2. r = -0.72 między temperaturą na zewnątrz a kosztami ogrzewania

   Odpowiedź: Silny ujemny związek. Wraz ze wzrostem temperatury koszty ogrzewania wyraźnie maleją.

3. r = 0.12 między rozmiarem buta a inteligencją

   Odpowiedź: Bardzo słaby/brak istotnego związku. Zmienne są praktycznie niezależne.

10.8 Najważniejsze rzeczy do zapamiętania (Important Points to Remember)

1. Korelacja mierzy siłę związku: Wartości od -1 do +1.

2. Korelacja ≠ przyczynowość (Correlation ≠ Causation): Wysoka korelacja nie dowodzi wpływu jednej zmiennej na drugą.

3. Dobierz właściwą metodę:

   • Pearson: Związki liniowe dla danych ciągłych.
   • Spearman: Związki monotoniczne lub dane rangowe.

4. Sprawdź założenia:

   • Pearson: liniowość i (w praktyce) rozkład zbliżony do normalnego.
   • Spearman: wymagana jedynie monotoniczność.

5. Uwaga na obserwacje odstające (outliers): Mogą silnie wpływać na korelację Pearsona.

6. Zawsze wizualizuj dane: Wykresy pomagają ocenić kształt zależności.

10.9 Podsumowanie: drzewko decyzyjne do analizy korelacji (Summary: Decision Tree for Correlation Analysis)

WYBÓR WŁAŚCIWEJ MIARY KORELACJI:

Czy dane są liczbowe (numeryczne)?
├─ TAK → Czy związek jest liniowy?
│   ├─ TAK → Użyj korelacji PEARSONA
│   └─ NIE → Czy związek jest monotoniczny?
│       ├─ TAK → Użyj korelacji SPEARMANA
│       └─ NIE → Rozważ metody nieliniowe
└─ NIE → Czy dane są porządkowe (rangi)?
    ├─ TAK → Użyj korelacji SPEARMANA
    └─ NIE → Użyj TABEL KRZYŻOWYCH dla danych kategorycznych

Ściąga (Quick Reference Card)

Miara	Kiedy używać (Use When)	Wzór (Formula)	Zakres (Range)
Kowariancja (Covariance)	Wstępne badanie związku	\displaystyle \frac{\sum_{i=1}^n (x_i-\bar{x})(y_i-\bar{y})}{n-1}	-\infty do +\infty
Pearson r	Związki liniowe, dane ciągłe	\displaystyle \frac{\operatorname{cov}(X,Y)}{s_X s_Y}	-1 do +1
Spearman \rho	Związki monotoniczne, rangi	\displaystyle 1-\frac{6\sum d_i^2}{n(n^2-1)}	-1 do +1
Tabele krzyżowe (Cross-tabs)	Zmienne kategoryczne	Zliczenia częstości	n/d

---

10.10 Analiza regresji OLS (Ordinary Least Squares): przewodnik na start (A Quick-start Guide)

Analiza regresji OLS: przewodnik na start

Wprowadzenie: czym jest analiza regresji?

Analiza regresji (regression analysis) pomaga zrozumieć i mierzyć zależności między obserwowalnymi wielkościami. To zestaw narzędzi matematycznych do identyfikowania wzorców w danych, które umożliwiają prognozowanie (prediction).

Rozważ pytania badawcze:

• Jak czas nauki wpływa na wynik testu?
• Jak doświadczenie wpływa na wynagrodzenie?
• Jak wydatki na reklamę oddziałują na sprzedaż?

Regresja dostarcza systematycznych metod, by na te pytania odpowiadać na podstawie realnych danych.

Punkt wyjścia: prosty przykład

Zacznijmy od konkretu. Zebrano dane o 20 studentach z Twojej klasy:

Student	Study Hours	Exam Score
Alex	2	68
Beth	4	74
Carlos	6	85
Diana	8	91
…	…	…

Po narysowaniu wykresu punktowego (scatter plot) chcesz znaleźć prostą, która najlepiej opisuje związek między godzinami nauki a wynikiem.

Ale co znaczy „najlepiej”? Właśnie to odkryjemy.

Dlaczego prawdziwe dane nie układają się w idealną linię

Zanim przejdziemy do rachunków, zrozummy, dlaczego punkty zwykle nie leżą na jednej prostej.

Modele deterministyczne vs. stochastyczne

Modele deterministyczne (deterministic models) opisują związki bez niepewności. Przykład z fizyki:

\text{Distance} = \text{Speed} \times \text{Time}

Jedziesz dokładnie 60 mph przez 2 godziny → zawsze 120 mil. Zero odchyleń.

Modele stochastyczne (stochastic models) uznają, że w danych naturalnie występuje losowość. Ogólna postać to:

Y = f(X) + \epsilon

Gdzie:

• Y — wielkość, którą prognozujemy (np. wynik testu),
• f(X) — wzorzec systematyczny (jak godziny nauki typowo wpływają na wyniki),
• \epsilon — „reszta”/szum: wszystko, czego nie mierzymy.

W naszym przykładzie dwoje studentów może uczyć się po 5 godzin, a jednak dostać różne oceny, bo:

• jedno lepiej spało,
• jedno ma talent do testów,
• jedno miało hałas na sali,
• pytania trafiły bardziej/mniej pod ich przygotowanie.

Ta losowość jest naturalna — tym zajmuje się \epsilon.

Prosty model regresji liniowej

Zależność między godzinami nauki a wynikiem zapisujemy jako:

Y_i = \beta_0 + \beta_1 X_i + \epsilon_i

Rozszyfrujmy:

• Y_i — wynik testu studenta i,
• X_i — godziny nauki studenta i,
• \beta_0 — wyraz wolny (intercept, „poziom bazowy” przy 0 godzin),
• \beta_1 — nachylenie (slope, przyrost punktów na godzinę),
• \epsilon_i — „wszystko inne” wpływające na wynik i.

Ważne: Prawdziwych wartości \beta_0, \beta_1 nie znamy. Szacujemy je z danych i oznaczamy „z daszkiem”: \hat{\beta}_0, \hat{\beta}_1.

Reszty: jak bardzo mylimy się w przewidywaniach?

Po dopasowaniu prostej możemy przewidzieć wyniki. Dla każdej obserwacji:

1. Wartość rzeczywista (y_i): faktyczny wynik,
2. Wartość przewidziana (\hat{y}_i): co „mówi” nasza prosta,
3. Reszta (e_i): różnica = Rzeczywista − Przewidziana.

Przykład:

Diana: 8 h nauki, wynik 91
Linia przewiduje: 88
Reszta: 91 − 88 = +3 (zaniżyliśmy)

Eric: 5 h nauki, wynik 70
Linia przewiduje: 79
Reszta: 70 − 79 = −9 (zawyżyliśmy)

Kluczowy pomysł: dlaczego reszty podnosimy do kwadratu?

Załóżmy reszty czterech studentów:

• A: +5
• B: −5
• C: +3
• D: −3

Suma: (+5) + (-5) + (+3) + (-3) = 0.

To nie znaczy, że przewidywania są idealne — błędy się znoszą.

Rozwiązanie: sumujemy kwadraty reszt:

• (+5)^2 = 25
• (-5)^2 = 25
• (+3)^2 = 9
• (-3)^2 = 9
• Suma kwadratów błędów = 68

Dlaczego to działa:

1. Brak znoszenia znaków, bo kwadraty są dodatnie,
2. Duże błędy ważą mocniej (10 punktów to 4× więcej niż 5 punktów),
3. Wygoda matematyczna: funkcje kwadratowe są gładkie i różniczkowalne.

Metoda zwykłych najmniejszych kwadratów (OLS)

OLS wybiera taką prostą, która minimalizuje sumę kwadratów reszt (SSE — Sum of Squared Errors):

\text{SSE} = \sum_{i=1}^n e_i^2 = \sum_{i=1}^n (y_i - \hat{y}_i)^2

Czyli:

\min_{\beta_0,\beta_1} \sum_{i=1}^n \big(y_i - (\beta_0 + \beta_1 x_i)\big)^2

„Po ludzku”: Znajdź takie \beta_0 i \beta_1, by łączny błąd (w kwadracie) przewidywań był jak najmniejszy.

Rozwiązanie matematyczne

Minimalizujemy SSE rachunkiem różniczkowym. Warunki pierwszego rzędu:

\frac{\partial SSE}{\partial \beta_0} = -2\sum_{i=1}^n (y_i - \beta_0 - \beta_1x_i) = 0

\frac{\partial SSE}{\partial \beta_1} = -2\sum_{i=1}^n x_i(y_i - \beta_0 - \beta_1x_i) = 0

Rozwiązanie układu:

\hat{\beta}_1 = \frac{\sum_{i=1}^n(x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n(x_i - \bar{x})^2} = \frac{\text{Cov}(X,Y)}{\text{Var}(X)}

\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}

Wnioski:

• Nachylenie zależy od tego, jak współzmieniają się X i Y (kowariancja) względem zmienności samego X (wariancja),
• Linia przechodzi przez punkt średnich (\bar{x}, \bar{y}).

Skąd wiemy, że linia jest „dobra”? Rozkład zmienności

Rozbijamy całkowitą zmienność wyników:

Całkowita suma kwadratów (SST — Total Sum of Squares)
„Jak bardzo ogólnie różnią się wyniki?”

SST = \sum_{i=1}^n (y_i - \bar{y})^2

Suma kwadratów regresji (SSR — Regression Sum of Squares)
„Ile zmienności wyjaśnia nasza linia?”

SSR = \sum_{i=1}^n (\hat{y}_i - \bar{y})^2

Suma kwadratów błędów (SSE — Error Sum of Squares)
„Ile nie wyjaśniamy?”

SSE = \sum_{i=1}^n (y_i - \hat{y}_i)^2

Tożsamość wariancyjna:

SST = SSR + SSE \quad\Rightarrow\quad \text{Całkowita} = \text{Wyjaśniona} + \text{Niewyjaśniona}

R^2: ocena dopasowania

Współczynnik determinacji (R^2) mówi, jaki odsetek zmienności wyjaśnia model:

R^2 = \frac{SSR}{SST} = 1 - \frac{SSE}{SST}

Interpretacja:

• R^2 = 0.75: „Godziny nauki wyjaśniają 75% zróżnicowania wyników”,
• R^2 = 0.30: „Model wyjaśnia 30% tego, czym różnią się wyniki”,
• R^2 = 1.00: perfekcyjne dopasowanie (w realu prawie nigdy),
• R^2 = 0.00: nie lepiej niż zgadywanie średniej.

Uwaga kontekstowa: w naukach społecznych 0.30 bywa świetne; w inżynierii oczekuje się bardzo wysokich R^2.

Interpretacja wyników

Załóżmy, że oszacowaliśmy \hat{\beta}_0 = 60 i \hat{\beta}_1 = 4.

Nachylenie (\hat{\beta}_1 = 4):

• „Każda dodatkowa godzina nauki wiąże się średnio z +4 punktami”,
• To efekt przeciętny, nie obietnica dla konkretnej osoby.

Wyraz wolny (\hat{\beta}_0 = 60):

• „Przy 0 godzinach nauki przewidujemy 60 punktów”,
• Często to tylko kotwica matematyczna — może nie mieć sensu praktycznego.

Równanie predykcji:

\widehat{\text{Wynik}} = 60 + 4 \times \text{Godziny nauki}

5 godzin → 60 + 4 \cdot 5 = 80 punktów.

Wielkość efektu i istotność praktyczna

Istotność statystyczna mówi, czy efekt istnieje; istotność praktyczna — czy ma znaczenie. Potrzebujemy obu.

Surowa wielkość efektu

To po prostu nachylenie \hat{\beta}_1.

Przykład: \hat{\beta}_1 = 4 pkt/godz.
Czy to „dużo”? Zależy od:

1. Skali wyniku (4/100 = 4% vs 4/500 = 0,8%),
2. Kosztu interwencji (czy 1h nauki warta 4 pkt?),
3. Progów decyzyjnych (czy 4 pkt zmienia ocenę?).

Standaryzowana wielkość efektu

Ułatwia porównania między badaniami:

\beta_{\text{std}} = \hat{\beta}_1 \times \frac{s_X}{s_Y}

gdzie s_X i s_Y to odchylenia standardowe X i Y.

Przykład: jeśli s_X = 2.5 h, s_Y = 12 pkt i \hat{\beta}_1 = 4:

\beta_{\text{std}} = 4 \cdot \frac{2.5}{12} = 0.83

„Przy wzroście X o 1 SD, Y rośnie o 0.83 SD”.

Wskazówki Cohena

Dla standaryzowanych współczynników:

• mały: |\beta| \approx 0.10 (~1% wariancji),
• średni: |\beta| \approx 0.30 (~9%),
• duży: |\beta| \approx 0.50 (~25%).

Dla R^2:

• mały: R^2 \approx 0.02,
• średni: R^2 \approx 0.13,
• duży: R^2 \approx 0.26.

Uwaga: to ogólne progi — normy różnią się między dziedzinami.

Przedziały ufności dla wielkości efektu

Dla surowego współczynnika:

CI = \hat{\beta}_1 \pm t_{\text{critical}} \cdot SE(\hat{\beta}_1)

Jeśli 95% CI = [3.2, 4.8], mówimy: „Z 95% pewnością prawdziwy efekt mieści się między 3.2 a 4.8 pkt/godz.”

Ocena istotności praktycznej

Weź pod uwagę:

1. Minimalnie istotną różnicę (MID),
2. Ile trzeba zmienić X, by osiągnąć sensowną zmianę Y,
3. Koszt-efektywność:

\text{Efektywność} = \frac{\text{Wielkość efektu}}{\text{Koszt interwencji}}

Przykład:
Efekt: 4 pkt/godz.
Próg zaliczenia: 70; średnia: 68 → 30 min nauki może zmienić niezal na zal → istotne praktycznie.

Niepewność

Szacunki pochodzą z próby, nie z całej populacji → niepewność.

Skąd niepewność?

• Masz 20 studentów, nie wszystkich,
• Próba może być nietypowa,
• Pomiary nie są doskonałe (raportowanie godzin nauki).

Przedziały ufności

Zamiast „efekt to dokładnie 4”, mówimy:

• „Szacujemy 4”,
• „95% CI: [3.2, 4.8]“.

Znaczenie: w wielu powtórzeniach 95% takich przedziałów zawiera prawdziwą wartość.

Testowanie istnienia zależności

Pytamy: „Gdyby prawdziwie nie było związku, jak mało prawdopodobny byłby obserwowany wzorzec?”

• p = 0.03: przy braku efektu tylko 3% szans na tak silny wzorzec „z przypadku”,
• p = 0.40: wzorzec „mógłby się łatwo zdarzyć” bez efektu.

Reguła kciuka: p < 0.05 → „statystycznie istotne”.

Gdy coś idzie źle: diagnostyka modelu

Szybkie wizualizacje:

1. Wykres punktowy: czy związek jest mniej więcej liniowy?
2. Reszty vs. przewidywania: powinien być losowy obłok,
3. Outliery: punkty bardzo odległe?

Sygnały ostrzegawcze:

• Wzorzec w resztach → brak liniowości lub zmienna pominięta,
• Rozszerzający się wachlarz reszt → heteroskedastyczność,
• Wpływowe obserwacje (influential) ciągną prostą,
• Pominięte zmienne → obciążenia (bias).

Założenia: kiedy OLS działa dobrze

1. Liniowość (linearity) — zależność w przybliżeniu prosta,
2. Niezależność (independence) — obserwacje od siebie niezależne,
3. Stała wariancja (homoskedastyczność) — rozrzut reszt podobny w całym zakresie,
4. Brak doskonałej współliniowości (no perfect multicollinearity) — w regresji wielorakiej predyktory nie są liniowo zależne „na 100%“,
5. Losowy dobór próby (random sampling) — dane reprezentatywne.

Podsumowanie

Co robi OLS:

• Dopasowuje prostą minimalizującą SSE,
• Szacuje, o ile średnio zmienia się Y przy zmianie X o 1,
• Podaje R^2 — ile zmienności wyjaśniamy,
• Kwantyfikuje niepewność (SE, CI, p-values).

Kroki praktyczne:

1. Wykres danych — czy linia ma sens?
2. Uruchom OLS → \hat{\beta}_0, \hat{\beta}_1,
3. Sprawdź R^2,
4. Oblicz wielkości efektu (surową i standaryzowaną),
5. Oceń istotność praktyczną,
6. Sprawdź przedziały ufności,
7. Obejrzyj reszty,
8. Decyduj, łącząc istotność statystyczną i praktyczną.

Kluczowe interpretacje / Key interpretations

Domyślny model: regresja OLS Y=\beta_0+\beta_1 X+\varepsilon (lub wieloraka: Y=\beta_0+\beta_1 X_1+\cdots+\beta_p X_p+\varepsilon).

• Nachylenie / Slope (\beta_1) PL: Przy wzroście X o 1 jednostkę (ceteris paribus), przeciętna wartość Y zmienia się o \beta_1 jednostek. ENG: When X increases by 1 unit (ceteris paribus), the expected value of Y changes by \beta_1 units.

• Standaryzowane nachylenie / Standardized slope \big(\beta_{1}^{(\mathrm{std})}\big) Definicja:

  \beta_{1}^{(\mathrm{std})} \;=\; \beta_1 \cdot \frac{s_X}{s_Y},

  gdzie s_X i s_Y to odchylenia standardowe X i Y. PL: Przy wzroście X o 1 odchylenie standardowe (SD), przeciętna wartość Y zmienia się o \beta_{1}^{(\mathrm{std})} odchyleń standardowych Y. ENG: For a 1 standard deviation (SD) increase in X, the expected value of Y changes by \beta_{1}^{(\mathrm{std})} SDs of Y. Uwaga/Note: W regresji prostej \beta_{1}^{(\mathrm{std})} = r (Pearson). / In simple regression, \beta_{1}^{(\mathrm{std})} = r (Pearson).

• Współczynnik determinacji / R^2 Definicja:

  R^2 \;=\; 1 - \frac{SS_{\mathrm{res}}}{SS_{\mathrm{tot}}}.

  PL: Model wyjaśnia 100\times R^2% zmienności Y względem modelu tylko z wyrazem wolnym (in-sample). ENG: The model explains 100\times R^2% of the variance in Y relative to the intercept-only model (in-sample). W wielu zmiennych rozważ: \text{adjusted } R^2. / With multiple predictors consider: adjusted R^2.

• Wartość p / P-value Formalnie/Formally:

  p \;=\; \Pr\!\big(\,|T|\ge |t_{\mathrm{obs}}| \mid H_0\,\big),

  gdzie T ma rozkład t przy H_0. PL: Zakładając prawdziwość H_0 i spełnione założenia modelu, prawdopodobieństwo uzyskania co najmniej tak ekstremalnej statystyki jak obserwowana wynosi p. ENG: Assuming H_0 and the model assumptions hold, p is the probability of observing a test statistic at least as extreme as the one obtained.

• Przedział ufności / Confidence interval (np. dla \beta_1) Konstrukcja/Construction:

  \hat{\beta}_1 \;\pm\; t_{1-\alpha/2,\ \mathrm{df}} \cdot \mathrm{SE}\!\left(\hat{\beta}_1\right).

  PL (ściśle): W długiej serii powtórzeń 95% tak skonstruowanych przedziałów zawiera prawdziwą wartość \beta_1; dla naszych danych oszacowanie mieści się w [\text{lower},\ \text{upper}]. ENG (strict): Over many repetitions, 95% of such intervals would contain the true \beta_1; for our data, the estimate lies within [\text{lower},\ \text{upper}]. PL (skrót dydaktyczny): „Jesteśmy 95% pewni, że \beta_1 leży w [\text{lower},\ \text{upper}].” ENG (teaching shorthand): “We are 95% confident that \beta_1 lies in [\text{lower},\ \text{upper}].”

Najczęstsze nieporozumienia / Common pitfalls

• PL: p nie jest prawdopodobieństwem, że H_0 jest prawdziwa. ENG: p is not the probability that H_0 is true.
• PL: 95% CI nie zawiera 95% obserwacji (od tego jest przedział predykcji). ENG: A 95% CI does not contain 95% of observations (that’s a prediction interval).
• PL/ENG: Wysokie R^2 ≠ przyczynowość / High R^2 ≠ causality. Zawsze sprawdzaj/Always check diagnozy reszt, skalę efektu, i dopasowanie poza próbą.

Pamiętaj:

• Asocjacja ≠ przyczynowość,
• Istotność statystyczna ≠ istotność praktyczna,
• „Każdy model jest błędny, niektóre są użyteczne”,
• Zawsze wizualizuj dane i reszty,
• Decyzje opieraj na wielkości efektu i niepewności.

OLS dostarcza uporządkowany, matematyczny sposób znajdowania wzorców w danych. Nie daje doskonałych prognoz, ale zapewnia najlepszą liniową aproksymację wraz z uczciwą oceną jej jakości i niepewności.

10.11 Ręczne obliczenia OLS krok po kroku

Badaczka chce zbadać zależność między godzinami nauki a wynikiem testu (6 studentów):

Student	Study Hours (X)	Exam Score (Y)
A	1	65
B	2	70
C	3	75
D	4	85
E	5	88
F	6	95

Celem jest wyznaczyć \hat{Y} = \hat{\beta}_0 + \hat{\beta}_1 X metodą OLS.

Krok 1: Średnie

Dla X:

\bar{X} = \frac{1+2+3+4+5+6}{6} = \frac{21}{6} = 3.5

Dla Y:

\bar{Y} = \frac{65+70+75+85+88+95}{6} = \frac{478}{6} = 79.67

Krok 2: Odchylenia od średnich

Student	X_i	Y_i	X_i - \bar{X}	Y_i - \bar{Y}
A	1	65	−2.5	−14.67
B	2	70	−1.5	−9.67
C	3	75	−0.5	−4.67
D	4	85	0.5	5.33
E	5	88	1.5	8.33
F	6	95	2.5	15.33

Krok 3: Iloczyny i kwadraty

Student	(X_i - \bar{X})(Y_i - \bar{Y})	(X_i - \bar{X})^2
A	36.68	6.25
B	14.51	2.25
C	2.34	0.25
D	2.67	0.25
E	12.50	2.25
F	38.33	6.25
Suma	107.03	17.50

Krok 4: Nachylenie \hat{\beta}_1

\hat{\beta}_1 = \frac{107.03}{17.50} = 6.12

Interpretacja: +6.12 punktu za każdą dodatkową godzinę nauki.

Krok 5: Wyraz wolny \hat{\beta}_0

\hat{\beta}_0 = 79.67 - 6.12 \cdot 3.5 = 58.25

Interpretacja: przy 0 godzinach przewidujemy 58.25.

Krok 6: Równanie regresji

\hat{Y} = 58.25 + 6.12 X

Krok 7: Predykcje i reszty

Student	X_i	Y_i	\hat{Y}_i = 58.25 + 6.12X_i	Reszta e_i = Y_i - \hat{Y}_i
A	1	65	64.37	0.63
B	2	70	70.49	−0.49
C	3	75	76.61	−1.61
D	4	85	82.73	2.27
E	5	88	88.85	−0.85
F	6	95	94.97	0.03

Kontrola: suma reszt ≈ 0 ✓

Krok 8: Sumy kwadratów

SST — całkowita:

Student	Y_i	(Y_i - \bar{Y})^2
A	65	215.21
B	70	93.51
C	75	21.81
D	85	28.41
E	88	69.39
F	95	235.01
Suma		SST = 663.34

SSR — wyjaśniona:

Student	\hat{Y}_i	(\hat{Y}_i - \bar{Y})^2
A	64.37	234.09
B	70.49	84.27
C	76.61	9.36
D	82.73	9.36
E	88.85	84.27
F	94.97	234.09
Suma		SSR = 655.44

SSE — błąd:

Student	e_i	e_i^2
A	0.63	0.40
B	−0.49	0.24
C	−1.61	2.59
D	2.27	5.15
E	−0.85	0.72
F	0.03	0.00
Suma		SSE = 9.10

Weryfikacja: SST \approx SSR + SSE
663.34 \approx 655.44 + 9.10 = 664.54 (drobne różnice zaokrągleń).

Krok 9: R^2

R^2 = \frac{SSR}{SST} = \frac{655.44}{663.34} = 0.988

Alternatywnie:

R^2 = 1 - \frac{SSE}{SST} = 1 - \frac{9.10}{663.34} = 0.986

Interpretacja: ~98.8% zmienności wyników wyjaśniają godziny nauki — bardzo silny związek.

Krok 10: Wielkości efektu

Surowa: 6.12 pkt/godz.

Standaryzowana:

• s_X = \sqrt{17.50/5} = 1.87,
• s_Y = \sqrt{663.34/5} = 11.52,
• \beta_{\text{std}} = 6.12 \cdot (1.87/11.52) = 0.99 → bardzo duży efekt (wg Cohena).

Krok 11: Istotność praktyczna

1. Skala: 6.12% na 100-punktowej skali / godz.,
2. Progi: zmiana oceny (10 pkt) ≈ 1.63 h,
3. Koszt-efekt: korzystny — sensowna inwestycja czasu.

Podsumowanie wyników

• Równanie: \hat{Y} = 58.25 + 6.12 X
• Nachylenie: 6.12 pkt/godz.
• Wyraz wolny: 58.25 pkt
• R^2: 0.988
• \beta_{\text{std}}: 0.99

W praktyce: każda godzina nauki to ≈ +6 pkt; dopasowanie znakomite; efekt istotny statystycznie i praktycznie.

Kontrola wyniku

Sprawdź, że linia przechodzi przez (\bar{X}, \bar{Y}):

58.25 + 6.12 \cdot 3.5 = 79.67 = \bar{Y} ✓

10.12 Kod R do weryfikacji obliczeń

# Krok 1: Dane
study_hours <- c(1, 2, 3, 4, 5, 6)      # X
exam_scores <- c(65, 70, 75, 85, 88, 95) # Y
n <- length(study_hours)

data <- data.frame(
  Student = LETTERS[1:6],
  X = study_hours,
  Y = exam_scores
)
print("Original Data:")

[1] "Original Data:"

print(data)

  Student X  Y
1       A 1 65
2       B 2 70
3       C 3 75
4       D 4 85
5       E 5 88
6       F 6 95

# Krok 2: Średnie
x_bar <- mean(study_hours)
y_bar <- mean(exam_scores)
cat("\nMeans:\n")

Means:

cat(sprintf("Mean of X (study hours): %.2f\n", x_bar))

Mean of X (study hours): 3.50

cat(sprintf("Mean of Y (exam scores): %.2f\n", y_bar))

Mean of Y (exam scores): 79.67

# Krok 3: Odchylenia od średnich
data$x_dev <- data$X - x_bar
data$y_dev <- data$Y - y_bar
cat("\nDeviations from means:\n")

Deviations from means:

print(data[, c("Student", "x_dev", "y_dev")])

  Student x_dev      y_dev
1       A  -2.5 -14.666667
2       B  -1.5  -9.666667
3       C  -0.5  -4.666667
4       D   0.5   5.333333
5       E   1.5   8.333333
6       F   2.5  15.333333

# Krok 4: Iloczyny i kwadraty
data$xy_product <- data$x_dev * data$y_dev
data$x_dev_sq <- data$x_dev^2

sum_xy_product <- sum(data$xy_product)
sum_x_dev_sq <- sum(data$x_dev_sq)
cat("\nSum of (Xi - X̄)(Yi - Ȳ):", round(sum_xy_product, 2), "\n")

Sum of (Xi - X̄)(Yi - Ȳ): 107 

cat("Sum of (Xi - X̄)²:", sum_x_dev_sq, "\n")

Sum of (Xi - X̄)²: 17.5 

# Krok 5: Nachylenie (beta_1)
beta_1_manual <- sum_xy_product / sum_x_dev_sq
cat("\nSlope (β₁) calculated manually:", round(beta_1_manual, 2), "\n")

Slope (β₁) calculated manually: 6.11 

# Krok 6: Wyraz wolny (beta_0)
beta_0_manual <- y_bar - beta_1_manual * x_bar
cat("Intercept (β₀) calculated manually:", round(beta_0_manual, 2), "\n")

Intercept (β₀) calculated manually: 58.27 

# Krok 7: Porównanie z lm()
model <- lm(Y ~ X, data = data)
cat("\n--- Verification with R's lm() function ---\n")

--- Verification with R's lm() function ---

print(summary(model)$coefficients[, 1:2])

             Estimate Std. Error
(Intercept) 58.266667  1.4045278
X            6.114286  0.3606495

# Krok 8: Predykcje i reszty
data$Y_hat <- beta_0_manual + beta_1_manual * data$X
data$residual <- data$Y - data$Y_hat

cat("\nPredicted values and residuals:\n")

Predicted values and residuals:

print(data[, c("Student", "X", "Y", "Y_hat", "residual")])

  Student X  Y    Y_hat    residual
1       A 1 65 64.38095  0.61904762
2       B 2 70 70.49524 -0.49523810
3       C 3 75 76.60952 -1.60952381
4       D 4 85 82.72381  2.27619048
5       E 5 88 88.83810 -0.83809524
6       F 6 95 94.95238  0.04761905

cat("Sum of residuals (should be ≈ 0):", round(sum(data$residual), 4), "\n")

Sum of residuals (should be ≈ 0): 0 

# Krok 9: Sumy kwadratów
SST <- sum((data$Y - y_bar)^2)
cat("\n--- Sum of Squares Calculations ---\n")

--- Sum of Squares Calculations ---

cat("SST (Total Sum of Squares):", round(SST, 2), "\n")

SST (Total Sum of Squares): 663.33 

SSR <- sum((data$Y_hat - y_bar)^2)
cat("SSR (Regression Sum of Squares):", round(SSR, 2), "\n")

SSR (Regression Sum of Squares): 654.23 

SSE <- sum(data$residual^2)
cat("SSE (Error Sum of Squares):", round(SSE, 2), "\n")

SSE (Error Sum of Squares): 9.1 

cat("\nVerification: SST = SSR + SSE\n")

Verification: SST = SSR + SSE

cat(sprintf("%.2f = %.2f + %.2f\n", SST, SSR, SSE))

663.33 = 654.23 + 9.10

cat("Difference (due to rounding):", round(SST - (SSR + SSE), 4), "\n")

Difference (due to rounding): 0 

# Krok 10: R-kwadrat
R_squared_method1 <- SSR / SST
R_squared_method2 <- 1 - (SSE / SST)
R_squared_lm <- summary(model)$r.squared

cat("\n--- R-squared Calculations ---\n")

--- R-squared Calculations ---

cat("R² (Method 1: SSR/SST):", round(R_squared_method1, 4), "\n")

R² (Method 1: SSR/SST): 0.9863 

cat("R² (Method 2: 1 - SSE/SST):", round(R_squared_method2, 4), "\n")

R² (Method 2: 1 - SSE/SST): 0.9863 

cat("R² (from lm function):", round(R_squared_lm, 4), "\n")

R² (from lm function): 0.9863 

# Krok 11: Wielkości efektu
cat("\n--- Effect Size Calculations ---\n")

--- Effect Size Calculations ---

cat("Raw effect size (slope):", round(beta_1_manual, 2), "points per hour\n")

Raw effect size (slope): 6.11 points per hour

sd_x <- sd(data$X)
sd_y <- sd(data$Y)
cat("\nStandard deviations:\n")

Standard deviations:

cat(sprintf("SD of X: %.2f\n", sd_x))

SD of X: 1.87

cat(sprintf("SD of Y: %.2f\n", sd_y))

SD of Y: 11.52

beta_std <- beta_1_manual * (sd_x / sd_y)
cat(sprintf("\nStandardized effect size: %.2f\n", beta_std))

Standardized effect size: 0.99

# Korelacja (dla regresji prostej |r| = sqrt(R^2))
correlation <- cor(data$X, data$Y)
cat("\nPearson correlation coefficient:", round(correlation, 4), "\n")

Pearson correlation coefficient: 0.9931 

cat("Square root of R²:", round(sqrt(R_squared_method1), 4), "\n")

Square root of R²: 0.9931 

cat("These should be equal (within rounding error)\n")

These should be equal (within rounding error)

# Krok 12: Wizualizacja
cat("\n--- Creating Visualization ---\n")

--- Creating Visualization ---

plot(data$X, data$Y, 
     main = "Study Hours vs Exam Scores with OLS Regression Line",
     xlab = "Study Hours",
     ylab = "Exam Score",
     pch = 19, col = "blue", cex = 1.5,
     xlim = c(0, 7), ylim = c(60, 100))

abline(a = beta_0_manual, b = beta_1_manual, col = "red", lwd = 2)

# Punkt średnich
points(x_bar, y_bar, pch = 15, col = "green", cex = 2)

# Reszty jako odcinki pionowe
for(i in 1:nrow(data)) {
  segments(data$X[i], data$Y[i], data$X[i], data$Y_hat[i], 
           col = "gray", lty = 2)
}

legend("bottomright", 
       legend = c("Observed Data", "Regression Line", "Mean Point", "Residuals"),
       col = c("blue", "red", "green", "gray"),
       pch = c(19, NA, 15, NA),
       lty = c(NA, 1, NA, 2),
       lwd = c(NA, 2, NA, 1))

# Równanie na wykresie
equation_text <- sprintf("Y = %.2f + %.2f*X", beta_0_manual, beta_1_manual)
r2_text <- sprintf("R² = %.3f", R_squared_method1)
text(1.5, 92, equation_text, pos = 4)
text(1.5, 89, r2_text, pos = 4)

Analiza regresji OLS - godziny nauki vs. wynik testu

# Podsumowanie końcowe
cat("\n========== FINAL SUMMARY ==========\n")

========== FINAL SUMMARY ==========

cat("Regression Equation: Y =", round(beta_0_manual, 2), "+", 
    round(beta_1_manual, 2), "* X\n")

Regression Equation: Y = 58.27 + 6.11 * X

cat("R-squared:", round(R_squared_method1, 4), 
    sprintf("(%.1f%% of variance explained)\n", R_squared_method1 * 100))

R-squared: 0.9863 (98.6% of variance explained)

cat("Standardized Effect Size:", round(beta_std, 2), "(Very Large Effect)\n")

Standardized Effect Size: 0.99 (Very Large Effect)

cat("\nInterpretation:\n")

Interpretation:

cat("- Each additional hour of study increases exam score by", 
    round(beta_1_manual, 2), "points\n")

- Each additional hour of study increases exam score by 6.11 points

cat("- Study hours explain", sprintf("%.1f%%", R_squared_method1 * 100), 
    "of the variation in exam scores\n")

- Study hours explain 98.6% of the variation in exam scores

cat("- The relationship is extremely strong (correlation =", 
    round(correlation, 3), ")\n")

- The relationship is extremely strong (correlation = 0.993 )

10.13 Jak uruchomić kod

1. Skopiuj cały blok kodu,
2. Wklej do RStudio,
3. Uruchom chunk po chunk lub cały dokument,
4. Porównaj wyniki z obliczeniami ręcznymi.

Co zobaczysz:

• Nachylenie: 6.12,
• Wyraz wolny: 58.25,
• R^2: ≈ 0.988,
• Efekt standaryzowany: ≈ 0.99,
• Wykres z punktami, linią regresji i resztami.

To potwierdza poprawność obliczeń manualnych.

---

10.14 Appendix A: Obliczanie Kowariancji, Korelacji Pearsona i Spearmana, oraz modelowanie OLS - przykład wprowadzający

Studentka politologii bada związek między wielkością okręgu wyborczego (DM) a wskaźnikiem dysproporcjonalności Gallaghera (GH) w wyborach parlamentarnych w 10 losowo wybranych demokracjach.

Dane dotyczące wielkości okręgu wyborczego (\text{DM}) i indeksu Gallaghera:

\text{DM} (X)	Gallagher (Y)
2	18,2
3	16,7
4	15,8
5	15,3
6	15,0
7	14,8
8	14,7
9	14,6
10	14,55
11	14,52

Krok 1: Obliczanie Podstawowych Statystyk

Obliczanie średnich:

Dla \text{DM} (X): \bar{X} = \frac{\sum_{i=1}^n X_i}{n}

Szczegółowe obliczenia:

2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 = 65 \bar{x} = \frac{65}{10} = 6,5

Dla indeksu Gallaghera (Y): \bar{Y} = \frac{\sum_{i=1}^n Y_i}{n}

Szczegółowe obliczenia:

18,2 + 16,7 + 15,8 + 15,3 + 15,0 + 14,8 + 14,7 + 14,6 + 14,55 + 14,52 = 154,17 \bar{y} = \frac{154,17}{10} = 15,417

Krok 2: Szczegółowe Obliczenia Kowariancji

Pełna tabela robocza ze wszystkimi obliczeniami:

i	X_i	Y_i	(X_i - \bar{X})	(Y_i - \bar{Y})	(X_i - \bar{X})(Y_i - \bar{Y})	(X_i - \bar{X})^2	(Y_i - \bar{Y})^2
1	2	18,2	-4,5	2,783	-12,5235	20,25	7,7451
2	3	16,7	-3,5	1,283	-4,4905	12,25	1,6461
3	4	15,8	-2,5	0,383	-0,9575	6,25	0,1467
4	5	15,3	-1,5	-0,117	0,1755	2,25	0,0137
5	6	15,0	-0,5	-0,417	0,2085	0,25	0,1739
6	7	14,8	0,5	-0,617	-0,3085	0,25	0,3807
7	8	14,7	1,5	-0,717	-1,0755	2,25	0,5141
8	9	14,6	2,5	-0,817	-2,0425	6,25	0,6675
9	10	14,55	3,5	-0,867	-3,0345	12,25	0,7517
10	11	14,52	4,5	-0,897	-4,0365	20,25	0,8047
Suma	65	154,17	0	0	-28,085	82,5	12,8442

Obliczanie kowariancji: \text{Cov}(X,Y) = \frac{\sum_{i=1}^n (X_i - \bar{X})(Y_i - \bar{Y})}{n-1}

\text{Cov}(X,Y) = \frac{-28,085}{9} = -3,120556

Krok 3: Obliczanie Odchylenia Standardowego

Dla \text{DM} (X): \sigma_X = \sqrt{\frac{\sum_{i=1}^n (X_i - \bar{X})^2}{n-1}}

\sigma_x = \sqrt{\frac{82,5}{9}} = \sqrt{9,1667} = 3,026582

Dla Gallaghera (Y): \sigma_Y = \sqrt{\frac{\sum_{i=1}^n (Y_i - \bar{Y})^2}{n-1}}

\sigma_y = \sqrt{\frac{12,8442}{9}} = \sqrt{1,4271} = 1,194612

Krok 4: Obliczanie Korelacji Pearsona

r = \frac{\text{Cov}(X,Y)}{\sigma_X \sigma_Y}

r = \frac{-3,120556}{3,026582 \times 1,194612} = \frac{-3,120556}{3,615752} = -0,863044

Krok 5: Obliczanie Korelacji Rangowej Spearmana

Pełna tabela rangowa ze wszystkimi obliczeniami:

i	X_i	Y_i	Ranga X_i	Ranga Y_i	d_i	d_i^2
1	2	18,2	1	10	-9	81
2	3	16,7	2	9	-7	49
3	4	15,8	3	8	-5	25
4	5	15,3	4	7	-3	9
5	6	15,0	5	6	-1	1
6	7	14,8	6	5	1	1
7	8	14,7	7	4	3	9
8	9	14,6	8	3	5	25
9	10	14,55	9	2	7	49
10	11	14,52	10	1	9	81
Suma						330

Obliczanie korelacji Spearmana: \rho = 1 - \frac{6\sum d_i^2}{n(n^2-1)}

\rho = 1 - \frac{6 \times 330}{10(100 - 1)} = 1 - \frac{1980}{990} = 1 - 2 = -1

Krok 6: Weryfikacja w R

# Tworzenie wektorów
DM <- c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11)
GH <- c(18.2, 16.7, 15.8, 15.3, 15.0, 14.8, 14.7, 14.6, 14.55, 14.52)

# Obliczanie kowariancji
cov(DM, GH)

[1] -3.120556

# Obliczanie korelacji
cor(DM, GH, method = "pearson")

[1] -0.8627742

cor(DM, GH, method = "spearman")

[1] -1

Krok 7: Podstawowa Wizualizacja

library(ggplot2)

# Tworzenie ramki danych
data <- data.frame(DM = DM, GH = GH)

# Tworzenie wykresu rozrzutu
ggplot(data, aes(x = DM, y = GH)) +
  geom_point(size = 3, color = "blue") +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(
    title = "Wielkość Okręgu vs Indeks Gallaghera",
    x = "Wielkość Okręgu (DM)",
    y = "Indeks Gallaghera (GH)"
  ) +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

Estymacja OLS i Miary Dopasowania Modelu

Krok 1: Obliczanie Estymatorów OLS

Korzystając z wcześniej obliczonych wartości:

• \sum(X_i - \bar{X})(Y_i - \bar{Y}) = -28,085
• \sum(X_i - \bar{X})^2 = 82,5
• \bar{X} = 6,5
• \bar{Y} = 15,417

Obliczanie nachylenia (\hat{\beta_1}):

\hat{\beta_1} = \frac{\sum(X_i - \bar{X})(Y_i - \bar{Y})}{\sum(X_i - \bar{X})^2}

\hat{\beta_1} = -28,085 ÷ 82,5 = -0,3404

Obliczanie wyrazu wolnego (\hat{\beta_0}): \hat{\beta_0} = \bar{Y} - \hat{\beta_1}\bar{X}

\hat{\beta_0} = 15,417 - (-0,3404 × 6,5) = 15,417 + 2,2126 = 17,6296

Zatem równanie regresji OLS ma postać: \hat{Y} = 17,6296 - 0,3404X

Krok 2: Obliczanie Wartości Dopasowanych i Reszt

Pełna tabela ze wszystkimi obliczeniami:

i	X_i	Y_i	\hat{Y}_i	e_i = Y_i - \hat{Y}_i	e_i^2	(Y_i - \bar{Y})^2	(\hat{Y}_i - \bar{Y})^2
1	2	18,2	16,9488	1,2512	1,5655	7,7451	2,3404
2	3	16,7	16,6084	0,0916	0,0084	1,6461	1,4241
3	4	15,8	16,2680	-0,4680	0,2190	0,1467	0,7225
4	5	15,3	15,9276	-0,6276	0,3939	0,0137	0,2601
5	6	15,0	15,5872	-0,5872	0,3448	0,1739	0,0289
6	7	14,8	15,2468	-0,4468	0,1996	0,3807	0,0290
7	8	14,7	14,9064	-0,2064	0,0426	0,5141	0,2610
8	9	14,6	14,5660	0,0340	0,0012	0,6675	0,7241
9	10	14,55	14,2256	0,3244	0,1052	0,7517	1,4184
10	11	14,52	13,8852	0,6348	0,4030	0,8047	2,3439
Suma	65	154,17	154,17	0	3,2832	12,8442	9,5524

Obliczenia dla wartości dopasowanych:

Dla X = 2:
Ŷ = 17,6296 + (-0,3404 × 2) = 16,9488

Dla X = 3:
Ŷ = 17,6296 + (-0,3404 × 3) = 16,6084

[... kontynuacja dla wszystkich wartości]

Krok 3: Obliczanie Miar Dopasowania

Suma kwadratów reszt (SSE): SSE = \sum e_i^2

SSE = 3,2832

Całkowita suma kwadratów (SST): SST = \sum(Y_i - \bar{Y})^2

SST = 12,8442

Suma kwadratów regresji (SSR): SSR = \sum(\hat{Y}_i - \bar{Y})^2

SSR = 9,5524

Weryfikacja dekompozycji: SST = SSR + SSE

12,8442 = 9,5524 + 3,2832 (w granicach błędu zaokrąglenia)

Obliczanie współczynnika determinacji R-kwadrat: R^2 = \frac{SSR}{SST} = 1 - \frac{SSE}{SST}

R² = 9,5524 ÷ 12,8442
   = 0,7438

Krok 4: Weryfikacja w R

# Dopasowanie modelu liniowego
model <- lm(GH ~ DM, data = data)

# Podsumowanie statystyk
summary(model)

Call:
lm(formula = GH ~ DM, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.62764 -0.46274 -0.08615  0.26624  1.25109 

Coefficients:
            Estimate Std. Error t value       Pr(>|t|)    
(Intercept) 17.62976    0.50121  35.174 0.000000000467 ***
DM          -0.34042    0.07053  -4.827        0.00131 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6406 on 8 degrees of freedom
Multiple R-squared:  0.7444,    Adjusted R-squared:  0.7124 
F-statistic:  23.3 on 1 and 8 DF,  p-value: 0.00131

# Ręczne obliczenie R-kwadrat
SST <- sum((GH - mean(GH))^2)
SSE <- sum(residuals(model)^2)
SSR <- SST - SSE
R2_manual <- SSR/SST
R2_manual

[1] 0.7443793

Krok 5: Analiza Reszt

# Tworzenie wykresów reszt
par(mfrow = c(2, 2))
plot(model)

Krok 6: Wykres Wartości Przewidywanych vs Rzeczywistych

# Tworzenie wykresu wartości przewidywanych vs rzeczywistych
ggplot(data.frame(
  Rzeczywiste = GH,
  Przewidywane = fitted(model)
), aes(x = Przewidywane, y = Rzeczywiste)) +
  geom_point(color = "blue", size = 3) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "red") +
  labs(
    title = "Wartości Przewidywane vs Rzeczywiste",
    x = "Przewidywany Indeks Gallaghera",
    y = "Rzeczywisty Indeks Gallaghera"
  ) +
  theme_minimal()

Modele z Transformacją Logarytmiczną

Krok 1: Transformacja Danych

Najpierw obliczamy logarytmy naturalne zmiennych:

i	X_i	Y_i	\ln(X_i)	\ln(Y_i)
1	2	18,2	0,6931	2,9014
2	3	16,7	1,0986	2,8154
3	4	15,8	1,3863	2,7600
4	5	15,3	1,6094	2,7278
5	6	15,0	1,7918	2,7081
6	7	14,8	1,9459	2,6946
7	8	14,7	2,0794	2,6878
8	9	14,6	2,1972	2,6810
9	10	14,55	2,3026	2,6777
10	11	14,52	2,3979	2,6757

Krok 2: Porównanie Różnych Specyfikacji Modelu

Szacujemy trzy alternatywne specyfikacje:

1. Model log-liniowy: \ln(Y_i) = \beta_0 + \beta_1 X_i + \epsilon_i
2. Model liniowo-logarytmiczny: Y_i = \beta_0 + \beta_1\ln(X_i) + \epsilon_i
3. Model log-log: \ln(Y_i) = \beta_0 + \beta_1\ln(X_i) + \epsilon_i

# Tworzenie zmiennych transformowanych
data$log_DM <- log(data$DM)
data$log_GH <- log(data$GH)

# Dopasowanie modeli
model_linear <- lm(GH ~ DM, data = data)
model_loglinear <- lm(log_GH ~ DM, data = data)
model_linearlog <- lm(GH ~ log_DM, data = data)
model_loglog <- lm(log_GH ~ log_DM, data = data)

# Porównanie wartości R-kwadrat
models_comparison <- data.frame(
  Model = c("Liniowy", "Log-liniowy", "Liniowo-logarytmiczny", "Log-log"),
  R_kwadrat = c(
    summary(model_linear)$r.squared,
    summary(model_loglinear)$r.squared,
    summary(model_linearlog)$r.squared,
    summary(model_loglog)$r.squared
  )
)

# Wyświetlenie porównania
models_comparison

                  Model R_kwadrat
1               Liniowy 0.7443793
2           Log-liniowy 0.7670346
3 Liniowo-logarytmiczny 0.9141560
4               Log-log 0.9288088

Krok 3: Porównanie Wizualne

# Tworzenie wykresów dla każdego modelu
p1 <- ggplot(data, aes(x = DM, y = GH)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Model Liniowy") +
  theme_minimal()

p2 <- ggplot(data, aes(x = DM, y = log_GH)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Model Log-liniowy") +
  theme_minimal()

p3 <- ggplot(data, aes(x = log_DM, y = GH)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Model Liniowo-logarytmiczny") +
  theme_minimal()

p4 <- ggplot(data, aes(x = log_DM, y = log_GH)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Model Log-log") +
  theme_minimal()

# Układanie wykresów w siatkę
library(gridExtra)
grid.arrange(p1, p2, p3, p4, ncol = 2)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Krok 4: Analiza Reszt dla Najlepszego Modelu

Na podstawie wartości R-kwadrat, analiza reszt dla najlepiej dopasowanego modelu:

# Wykresy reszt dla najlepszego modelu
par(mfrow = c(2, 2))
plot(model_linearlog)

Krok 5: Interpretacja Najlepszego Modelu

Współczynniki modelu liniowo-logarytmicznego:

summary(model_linearlog)

Call:
lm(formula = GH ~ log_DM, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.40702 -0.30207 -0.04907  0.22905  0.60549 

Coefficients:
            Estimate Std. Error t value        Pr(>|t|)    
(Intercept)  19.0223     0.4079   46.64 0.0000000000494 ***
log_DM       -2.0599     0.2232   -9.23 0.0000153880425 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3712 on 8 degrees of freedom
Multiple R-squared:  0.9142,    Adjusted R-squared:  0.9034 
F-statistic: 85.19 on 1 and 8 DF,  p-value: 0.00001539

Interpretacja:

• \hat{\beta_0} reprezentuje oczekiwany Indeks Gallaghera, gdy ln(DM) = 0 (czyli gdy DM = 1)
• \hat{\beta_1} reprezentuje zmianę Indeksu Gallaghera związaną z jednostkowym wzrostem ln(DM)

Krok 6: Predykcje Modelu

# Tworzenie wykresu predykcji dla najlepszego modelu
ggplot(data, aes(x = log_DM, y = GH)) +
  geom_point(color = "blue", size = 3) +
  geom_smooth(method = "lm", se = TRUE, color = "red") +
  labs(
    title = "Model Liniowo-logarytmiczny: Indeks Gallaghera vs ln(Wielkość Okręgu)",
    x = "ln(Wielkość Okręgu)",
    y = "Indeks Gallaghera"
  ) +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

Krok 7: Analiza Elastyczności

Dla modelu log-log współczynniki bezpośrednio reprezentują elastyczności. Obliczenie średniej elastyczności dla modelu liniowo-logarytmicznego:

# Obliczenie elastyczności przy wartościach średnich
mean_DM <- mean(data$DM)
mean_GH <- mean(data$GH)
beta1 <- coef(model_linearlog)[2]
elastycznosc <- beta1 * (1/mean_GH)
elastycznosc

    log_DM 
-0.1336136 

Wartość ta reprezentuje procentową zmianę Indeksu Gallaghera przy jednoprocentowej zmianie Wielkości Okręgu.

---

10.15 Appendix B. Prosty Model Regresji Liniowej (MNK/OLS) - przykłady różne

10.16 Przykład 0. Czas Nauki vs. Ocena: Przewodnik po MNK i R-kwadrat

library(ggplot2)
library(dplyr)
library(gridExtra)

# Nasz prosty zbiór danych
practice <- c(1, 2, 3, 4, 5, 6)
skill <- c(3, 7, 5, 8, 10, 9)

# Tworzenie ramki danych
data <- data.frame(
  Student = 1:6,
  Practice = practice,
  Skill = skill
)

print(data)

  Student Practice Skill
1       1        1     3
2       2        2     7
3       3        3     5
4       4        4     8
5       5        5    10
6       6        6     9

Obliczenia Ręczne

Krok 1: Obliczenie Średnich

Średnia to wartość przeciętna, obliczana przez zsumowanie wszystkich obserwacji i podzielenie przez ich liczbę.

Średnia Godzin Ćwiczeń (X̄):

\bar{X} = \frac{\sum_{i=1}^{n} X_i}{n} = \frac{1 + 2 + 3 + 4 + 5 + 6}{6} = \frac{21}{6} = 3.5

Średnia Ocen Umiejętności (Ȳ):

\bar{Y} = \frac{\sum_{i=1}^{n} Y_i}{n} = \frac{3 + 7 + 5 + 8 + 10 + 9}{6} = \frac{42}{6} = 7

# Obliczenie średnich
mean_x <- mean(practice)
mean_y <- mean(skill)

cat("Średnia Godzin Ćwiczeń:", mean_x, "\n")

Średnia Godzin Ćwiczeń: 3.5 

cat("Średnia Oceny Umiejętności:", mean_y, "\n")

Średnia Oceny Umiejętności: 7 

Krok 2: Obliczenie Wariancji i Odchylenia Standardowego

Wariancja mierzy, jak bardzo dane są rozproszone wokół średniej. Używamy wzoru na wariancję próbkową (dzielenie przez n-1).

Wariancja X (Godziny Ćwiczeń):

Najpierw obliczamy odchylenia od średniej (X_i - \bar{X}):

Student	X_i	X_i - \bar{X}	(X_i - \bar{X})^2
1	1	1 - 3.5 = -2.5	(-2.5)^2 = 6.25
2	2	2 - 3.5 = -1.5	(-1.5)^2 = 2.25
3	3	3 - 3.5 = -0.5	(-0.5)^2 = 0.25
4	4	4 - 3.5 = 0.5	(0.5)^2 = 0.25
5	5	5 - 3.5 = 1.5	(1.5)^2 = 2.25
6	6	6 - 3.5 = 2.5	(2.5)^2 = 6.25
Suma			17.5

s^2_X = \frac{\sum(X_i - \bar{X})^2}{n-1} = \frac{17.5}{5} = 3.5

s_X = \sqrt{3.5} = 1.871

Wariancja Y (Oceny Umiejętności):

Student	Y_i	Y_i - \bar{Y}	(Y_i - \bar{Y})^2
1	3	3 - 7 = -4	(-4)^2 = 16
2	7	7 - 7 = 0	(0)^2 = 0
3	5	5 - 7 = -2	(-2)^2 = 4
4	8	8 - 7 = 1	(1)^2 = 1
5	10	10 - 7 = 3	(3)^2 = 9
6	9	9 - 7 = 2	(2)^2 = 4
Suma			34

s^2_Y = \frac{\sum(Y_i - \bar{Y})^2}{n-1} = \frac{34}{5} = 6.8

s_Y = \sqrt{6.8} = 2.608

# Weryfikacja wariancji i odchylenia standardowego
cat("Wariancja Ćwiczeń:", var(practice), "\n")

Wariancja Ćwiczeń: 3.5 

cat("Odch. Stand. Ćwiczeń:", sd(practice), "\n")

Odch. Stand. Ćwiczeń: 1.870829 

cat("Wariancja Umiejętności:", var(skill), "\n")

Wariancja Umiejętności: 6.8 

cat("Odch. Stand. Umiejętności:", sd(skill), "\n")

Odch. Stand. Umiejętności: 2.607681 

Krok 3: Obliczenie Kowariancji

Kowariancja mierzy, jak dwie zmienne zmieniają się wspólnie. Dodatnia kowariancja wskazuje, że gdy jedna zmienna rośnie, druga również ma tendencję do wzrostu.

s_{XY} = \frac{\sum_{i=1}^{n}(X_i - \bar{X})(Y_i - \bar{Y})}{n-1}

Obliczmy iloczyny dla każdej obserwacji:

Student	(X_i - \bar{X})	(Y_i - \bar{Y})	(X_i - \bar{X})(Y_i - \bar{Y})
1	-2.5	-4	(-2.5) \times (-4) = 10.0
2	-1.5	0	(-1.5) \times (0) = 0.0
3	-0.5	-2	(-0.5) \times (-2) = 1.0
4	0.5	1	(0.5) \times (1) = 0.5
5	1.5	3	(1.5) \times (3) = 4.5
6	2.5	2	(2.5) \times (2) = 5.0
Suma			21.0

s_{XY} = \frac{21.0}{5} = 4.2

# Weryfikacja kowariancji
cat("Kowariancja:", cov(practice, skill), "\n")

Kowariancja: 4.2 

Krok 4: Obliczenie Współczynnika Korelacji Pearsona

Współczynnik korelacji standaryzuje kowariancję do skali od -1 do +1:

r = \frac{s_{XY}}{s_X \cdot s_Y} = \frac{4.2}{1.871 \times 2.608} = \frac{4.2}{4.879} = 0.861

Otrzymujemy korelację wynoszącą 0.861, co wskazuje na silny dodatni związek między godzinami ćwiczeń a oceną umiejętności.

# Weryfikacja korelacji
cat("Korelacja:", cor(practice, skill), "\n")

Korelacja: 0.8609161 

Krok 5: Obliczenie Współczynników Regresji MNK

Metoda Najmniejszych Kwadratów (MNK) znajduje wartości \beta_0 i \beta_1, które minimalizują sumę kwadratów błędów.

Estymator nachylenia:

\hat{\beta_1} = \frac{s_{XY}}{s^2_X} = \frac{\sum(X_i - \bar{X})(Y_i - \bar{Y})}{\sum(X_i - \bar{X})^2}

Używając naszych obliczonych wartości:

\hat{\beta_1} = \frac{4.2}{3.5} = 1.2

Estymator wyrazu wolnego:

\hat{\beta_0} = \bar{Y} - \hat{\beta_1}\bar{X}

\hat{\beta_0} = 7 - (1.2 \times 3.5) = 7 - 4.2 = 2.8

Nasze równanie regresji:

\hat{Y} = 2.8 + 1.2X

To oznacza:

• Gdy godziny ćwiczeń = 0, przewidywana umiejętność = 2.8

• Każda dodatkowa godzina ćwiczeń zwiększa ocenę umiejętności o 1.2 punktu

# Dopasowanie modelu
model <- lm(skill ~ practice)
coef_model <- coef(model)

cat("Wyraz wolny (β₀):", coef_model[1], "\n")

Wyraz wolny (β₀): 2.8 

cat("Nachylenie (β₁):", coef_model[2], "\n")

Nachylenie (β₁): 1.2 

# Obliczenie przewidywań
data$predicted <- predict(model)
data$residual <- residuals(model)

Krok 6: Obliczenie Wartości Przewidywanych i Reszt

Używając \hat{Y} = 2.8 + 1.2X:

Student	X_i	Y_i	\hat{Y}_i = 2.8 + 1.2X_i	Reszta (Y_i - \hat{Y}_i)
1	1	3	2.8 + 1.2(1) = 4.0	3 - 4.0 = -1.0
2	2	7	2.8 + 1.2(2) = 5.2	7 - 5.2 = 1.8
3	3	5	2.8 + 1.2(3) = 6.4	5 - 6.4 = -1.4
4	4	8	2.8 + 1.2(4) = 7.6	8 - 7.6 = 0.4
5	5	10	2.8 + 1.2(5) = 8.8	10 - 8.8 = 1.2
6	6	9	2.8 + 1.2(6) = 10.0	9 - 10.0 = -1.0

Krok 7: Obliczenie Sum Kwadratów

SST (Całkowita Suma Kwadratów) - Całkowita zmienność Y:

SST = \sum(Y_i - \bar{Y})^2

Z naszych wcześniejszych obliczeń wariancji:

SST = (n-1) \times s^2_Y = 5 \times 6.8 = 34

SSE (Suma Kwadratów Błędów) - Zmienność niewyjaśniona:

SSE = \sum(Y_i - \hat{Y}_i)^2

Student	Y_i	\hat{Y}_i	(Y_i - \hat{Y}_i)	(Y_i - \hat{Y}_i)^2
1	3	4.0	-1.0	1.00
2	7	5.2	1.8	3.24
3	5	6.4	-1.4	1.96
4	8	7.6	0.4	0.16
5	10	8.8	1.2	1.44
6	9	10.0	-1.0	1.00
Suma				8.80

SSE = 8.80

SSR (Suma Kwadratów Regresji) - Zmienność wyjaśniona:

SSR = SST - SSE = 34 - 8.80 = 25.20

Krok 8: Obliczenie R-kwadrat

R-kwadrat informuje nas, jaka część całkowitej zmienności Y jest wyjaśniona przez nasz model:

R^2 = \frac{SSR}{SST} = \frac{25.20}{34} = 0.741

Alternatywny wzór:

R^2 = 1 - \frac{SSE}{SST} = 1 - \frac{8.80}{34} = 1 - 0.259 = 0.741

Lub po prostu:

R^2 = r^2 = (0.861)^2 = 0.741

To oznacza, że nasz model wyjaśnia 74.1% zmienności ocen umiejętności!

# Weryfikacja sum kwadratów i R-kwadrat
SST <- sum((skill - mean_y)^2)
SSR <- sum((data$predicted - mean_y)^2)
SSE <- sum(data$residual^2)
r_squared <- SSR / SST

cat("SST (Całkowita):", SST, "\n")

SST (Całkowita): 34 

cat("SSR (Wyjaśniona):", SSR, "\n")

SSR (Wyjaśniona): 25.2 

cat("SSE (Niewyjaśniona):", SSE, "\n")

SSE (Niewyjaśniona): 8.8 

cat("R-kwadrat:", r_squared, "\n")

R-kwadrat: 0.7411765 

cat("R-kwadrat (z korelacji):", cor(practice, skill)^2, "\n")

R-kwadrat (z korelacji): 0.7411765 

Wizualizacja 1: Linia Najlepszego Dopasowania MNK

Ten wykres pokazuje, jak MNK minimalizuje sumę kwadratów reszt (pionowe odległości od punktów do linii).

ggplot(data, aes(x = Practice, y = Skill)) +
  geom_hline(yintercept = mean_y, linetype = "dashed", 
             color = "gray50", linewidth = 0.8, alpha = 0.7) +
  geom_segment(aes(xend = Practice, yend = predicted), 
               color = "red", linetype = "dotted", linewidth = 1) +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linewidth = 1.2) +
  geom_point(size = 5, color = "darkblue") +
  geom_point(aes(y = predicted), size = 3, color = "blue", shape = 17) +
  annotate("text", x = 1.5, y = mean_y + 0.3, 
           label = paste0("Średnia (ȳ = ", mean_y, ")"), 
           color = "gray30", size = 4) +
  annotate("text", x = 5, y = 4, 
           label = "Reszty (błędy)\nMNK minimalizuje Σ(reszty²)", 
           color = "red", size = 3.5, hjust = 0) +
  annotate("text", x = 2, y = 10.5, 
           label = paste0("ŷ = ", round(coef_model[1], 1), 
                         " + ", round(coef_model[2], 1), "x"), 
           color = "blue", size = 5, fontface = "bold") +
  labs(
    title = "Regresja MNK: Minimalizacja Kwadratów Reszt",
    subtitle = "Niebieskie trójkąty to wartości przewidywane; czerwone linie pokazują reszty",
    x = "Godziny Ćwiczeń",
    y = "Ocena Umiejętności"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", size = 14))

`geom_smooth()` using formula = 'y ~ x'

Kluczowa Obserwacja: MNK znajduje unikalną linię, która czyni sumę kwadratów czerwonych odległości jak najmniejszą!

Wizualizacja 2: Dekompozycja Wariancji (SST = SSR + SSE)

To pokazuje, jak całkowita zmienność jest dzielona na składniki wyjaśniony i niewyjaśniony.

# Tworzenie wykresu dekompozycji
ggplot(data, aes(x = Practice)) +
  # Całkowite odchylenie (SST)
  geom_segment(aes(y = mean_y, yend = Skill, xend = Practice), 
               color = "purple", linewidth = 1.2, alpha = 0.6,
               arrow = arrow(length = unit(0.15, "inches"))) +
  
  # Odchylenie wyjaśnione (SSR)
  geom_segment(aes(y = mean_y, yend = predicted, xend = Practice), 
               color = "green", linewidth = 1.5, alpha = 0.8) +
  
  # Odchylenie niewyjaśnione (SSE)
  geom_segment(aes(y = predicted, yend = Skill, xend = Practice), 
               color = "red", linewidth = 1.2, alpha = 0.8,
               linetype = "dashed") +
  
  # Linia średniej
  geom_hline(yintercept = mean_y, linetype = "solid", 
             color = "gray40", linewidth = 1) +
  
  # Linia regresji
  geom_smooth(aes(y = Skill), method = "lm", se = FALSE, 
              color = "blue", linewidth = 1) +
  
  # Punkty
  geom_point(aes(y = Skill), size = 5, color = "darkblue") +
  geom_point(aes(y = predicted), size = 3, color = "blue", shape = 15) +
  
  # Adnotacje
  annotate("text", x = 6.5, y = mean_y, 
           label = "Średnia", color = "gray40", size = 4, hjust = 0) +
  annotate("text", x = 6.5, y = 9.5, 
           label = "Linia Regresji", color = "blue", size = 4, hjust = 0) +
  
  # Legenda
  annotate("segment", x = 0.5, xend = 0.5, y = 2, yend = 3.5,
           color = "purple", linewidth = 1.2, 
           arrow = arrow(length = unit(0.12, "inches"))) +
  annotate("text", x = 0.7, y = 2.75, 
           label = "Całkowite (SST)", color = "purple", 
           size = 3.5, hjust = 0, fontface = "bold") +
  
  annotate("segment", x = 0.5, xend = 0.5, y = 4.5, yend = 5.5,
           color = "green", linewidth = 1.5) +
  annotate("text", x = 0.7, y = 5, 
           label = "Wyjaśnione (SSR)", color = "green", 
           size = 3.5, hjust = 0, fontface = "bold") +
  
  annotate("segment", x = 0.5, xend = 0.5, y = 6.5, yend = 7.5,
           color = "red", linewidth = 1.2, linetype = "dashed") +
  annotate("text", x = 0.7, y = 7, 
           label = "Niewyjaśnione (SSE)", color = "red", 
           size = 3.5, hjust = 0, fontface = "bold") +
  
  labs(
    title = "Dekompozycja Wariancji: SST = SSR + SSE",
    subtitle = "Fioletowe = odchylenie całkowite | Zielone = wyjaśnione przez model | Czerwone = błąd resztowy",
    x = "Godziny Ćwiczeń",
    y = "Ocena Umiejętności"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", size = 14)) +
  coord_cartesian(xlim = c(0.3, 7.5), ylim = c(2, 11))

`geom_smooth()` using formula = 'y ~ x'

Dekompozycja Matematyczna:

Dla każdej obserwacji: (Y_i - \bar{Y})^2 = (\hat{Y}_i - \bar{Y})^2 + (Y_i - \hat{Y}_i)^2

• Fioletowe strzałki: Całkowite odchylenie od średniej = (Y_i - \bar{Y})

• Zielone słupki: Wyjaśnienie modelu = (\hat{Y}_i - \bar{Y})

• Czerwone przerywane: Co zostaje = (Y_i - \hat{Y}_i)

Wizualizacja 3: R-kwadrat jako Proporcja

# Obliczenie sum kwadratów
SST <- sum((skill - mean_y)^2)
SSR <- sum((data$predicted - mean_y)^2)
SSE <- sum(data$residual^2)
r_squared <- SSR / SST

# Dane do wykresu słupkowego
ss_data <- data.frame(
  Component = c("Całkowita (SST)", "Wyjaśniona (SSR)", "Niewyjaśniona (SSE)"),
  Value = c(SST, SSR, SSE),
  Color = c("purple", "green", "red")
)

p1 <- ggplot(ss_data, aes(x = Component, y = Value, fill = Component)) +
  geom_bar(stat = "identity", alpha = 0.7, color = "black") +
  geom_text(aes(label = round(Value, 2)), vjust = -0.5, size = 5, fontface = "bold") +
  scale_fill_manual(values = c("green", "red", "purple")) +
  labs(
    title = "Rozkład Sum Kwadratów",
    y = "Suma Kwadratów",
    x = ""
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold"))

# Wizualizacja proporcji
prop_data <- data.frame(
  Component = c("Wyjaśniona\n(SSR)", "Niewyjaśniona\n(SSE)"),
  Proportion = c(SSR/SST, SSE/SST),
  Percentage = c(SSR/SST * 100, SSE/SST * 100)
)

p2 <- ggplot(prop_data, aes(x = "", y = Proportion, fill = Component)) +
  geom_bar(stat = "identity", width = 1, color = "white", linewidth = 2) +
  coord_polar(theta = "y") +
  scale_fill_manual(values = c("green3", "red3")) +
  geom_text(aes(label = paste0(round(Percentage, 1), "%")), 
            position = position_stack(vjust = 0.5), 
            size = 6, fontface = "bold", color = "white") +
  labs(title = paste0("R² = ", round(r_squared, 3))) +
  theme_void(base_size = 12) +
  theme(plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
        legend.position = "bottom",
        legend.title = element_blank())

grid.arrange(p1, p2, ncol = 2)

Wzory na R-kwadrat:

R^2 = \frac{SSR}{SST} = \frac{25.20}{34} = 0.74

R^2 = 1 - \frac{SSE}{SST} = 1 - \frac{8.80}{34} = 0.74

Interpretacja: Nasz model wyjaśnia 74% zmienności ocen umiejętności!

Wizualizacja 4: R² jako Korelacja między Odchyleniami

To pokazuje geometryczną interpretację R²: jak dobrze przewidywane odchylenia od średniej pasują do rzeczywistych odchyleń.

# Obliczenie odchyleń
data$actual_dev <- skill - mean_y
data$predicted_dev <- data$predicted - mean_y

# Porównanie obok siebie
p1 <- ggplot(data, aes(x = Practice)) +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.8) +
  geom_segment(aes(y = 0, yend = actual_dev, xend = Practice), 
               color = "purple", linewidth = 2, alpha = 0.7,
               arrow = arrow(length = unit(0.15, "inches"))) +
  geom_point(aes(y = actual_dev), size = 4, color = "purple") +
  labs(
    title = "Rzeczywiste Odchylenia od Średniej",
    subtitle = expression(Y[i] - bar(Y)),
    x = "Godziny Ćwiczeń",
    y = "Odchylenie od Średniej"
  ) +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

p2 <- ggplot(data, aes(x = Practice)) +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.8) +
  geom_segment(aes(y = 0, yend = predicted_dev, xend = Practice), 
               color = "green", linewidth = 2, alpha = 0.7,
               arrow = arrow(length = unit(0.15, "inches"))) +
  geom_point(aes(y = predicted_dev), size = 4, color = "green") +
  labs(
    title = "Przewidywane Odchylenia od Średniej",
    subtitle = expression(hat(Y)[i] - bar(Y)),
    x = "Godziny Ćwiczeń",
    y = "Odchylenie od Średniej"
  ) +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold"))

grid.arrange(p1, p2, ncol = 2)

Wizualizacja 5: Wykres Rozrzutu Odchyleń (Interpretacja R²)

ggplot(data, aes(x = predicted_dev, y = actual_dev)) +
  geom_abline(slope = 1, intercept = 0, 
              linetype = "dashed", color = "gray50", linewidth = 1) +
  geom_point(size = 5, color = "darkblue", alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE, color = "blue", 
              fill = "lightblue", alpha = 0.3) +
  annotate("text", x = -3, y = 2.5, 
           label = paste0("Gdyby dopasowanie idealne:\nwszystkie punkty na tej linii\n(R² = 1)"), 
           color = "gray40", size = 3.5, hjust = 0) +
  annotate("text", x = 1.5, y = -3, 
           label = paste0("Korelacja = ", round(sqrt(r_squared), 3), 
                         "\nR² = ", round(r_squared, 3)), 
           color = "blue", size = 5, fontface = "bold") +
  labs(
    title = "R² Mierzy, Jak Dobrze Przewidywane Odchylenia Pasują do Rzeczywistych",
    subtitle = "Każdy punkt porównuje przewidywanie modelu z rzeczywistością (obie względem średniej)",
    x = expression(paste("Odchylenie Przewidywane: ", hat(Y)[i] - bar(Y))),
    y = expression(paste("Odchylenie Rzeczywiste: ", Y[i] - bar(Y)))
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", size = 13))

`geom_smooth()` using formula = 'y ~ x'

Kluczowa Obserwacja:

R² to dosłownie kwadrat korelacji między:

• Tym, co przewiduje model (względem średniej)

• Tym, co rzeczywiście się zdarzyło (względem średniej)

Jeśli R^2 = 1: wszystkie punkty leżą dokładnie na przekątnej (idealne przewidywania)

Jeśli R^2 = 0: brak związku między przewidywanymi a rzeczywistymi odchyleniami

Tabela Podsumowująca

# Tworzenie tabeli podsumowującej
summary_stats <- data.frame(
  Statystyka = c("Liczebność próby (n)", 
                "Średnia Ćwiczeń (X̄)", 
                "Średnia Umiejętności (Ȳ)",
                "Korelacja (r)",
                "Wyraz wolny (β₀)",
                "Nachylenie (β₁)",
                "R-kwadrat",
                "SST (Całkowita)",
                "SSR (Wyjaśniona)",
                "SSE (Niewyjaśniona)"),
  Wartość = c(6,
            round(mean_x, 2),
            round(mean_y, 2),
            round(cor(practice, skill), 3),
            round(coef_model[1], 2),
            round(coef_model[2], 2),
            round(r_squared, 3),
            round(SST, 2),
            round(SSR, 2),
            round(SSE, 2))
)

knitr::kable(summary_stats, align = c("l", "r"))

Statystyka	Wartość
Liczebność próby (n)	6.000
Średnia Ćwiczeń (X̄)	3.500
Średnia Umiejętności (Ȳ)	7.000
Korelacja (r)	0.861
Wyraz wolny (β₀)	2.800
Nachylenie (β₁)	1.200
R-kwadrat	0.741
SST (Całkowita)	34.000
SSR (Wyjaśniona)	25.200
SSE (Niewyjaśniona)	8.800

Kluczowe Wnioski

Regresja MNK:

• Znajduje linię minimalizującą sumę kwadratów reszt

• Produkuje nieobciążone estymatory o minimalnej wariancji (BLUE)

R-kwadrat (0.74) oznacza:

1. 74% zmienności umiejętności jest wyjaśnione przez godziny ćwiczeń

2. Korelacja między przewidywanymi a rzeczywistymi odchyleniami wynosi \sqrt{0.74} = 0.86

3. SSR stanowi 74% SST; SSE stanowi 26% SST

Interpretacja Geometryczna:

• Całkowita zmienność = odległość każdego punktu od średniej

• Model uchwytuje 74% tych odległości przez linię regresji

• Pozostałe 26% jest niewyjaśnione (reszty)

Implikacja Praktyczna:

Każda dodatkowa godzina ćwiczeń zwiększa oczekiwaną umiejętność o 1.2 punktu, a ten związek wyjaśnia większość (ale nie całość) obserwowanej zmienności!

10.17 Przykład 1. Analiza związku między dobrobytem ekonomicznym a frekwencją wyborczą

Analiza związku między dobrobytem ekonomicznym a frekwencją wyborczą w dzielnicach średniej wielkości europejskiego miasta na podstawie danych z wyborów samorządowych.

Dane

Próba obejmuje pięć reprezentatywnych dzielnic:

Dzielnica	Dochód (tys. €)	Frekwencja (%)
Dzielnica A	50	60
Dzielnica B	45	56
Dzielnica C	56	70
Dzielnica D	40	50
Dzielnica E	60	75

# Wczytanie bibliotek
library(tidyverse)

# Utworzenie zbioru danych
dane <- data.frame(
  dzielnica = c("Dzielnica A", "Dzielnica B", "Dzielnica C", "Dzielnica D", "Dzielnica E"),
  dochod = c(50, 45, 56, 40, 60),
  frekwencja = c(60, 56, 70, 50, 75)
)

# Podgląd danych
dane

    dzielnica dochod frekwencja
1 Dzielnica A     50         60
2 Dzielnica B     45         56
3 Dzielnica C     56         70
4 Dzielnica D     40         50
5 Dzielnica E     60         75

Część 1: Statystyki opisowe

# Statystyki dla dochodu
cat("DOCHÓD (tys. €):\n")

DOCHÓD (tys. €):

cat("Średnia:", mean(dane$dochod), "\n")

Średnia: 50.2 

cat("Mediana:", median(dane$dochod), "\n")

Mediana: 50 

cat("Odchylenie standardowe:", round(sd(dane$dochod), 2), "\n")

Odchylenie standardowe: 8.07 

cat("Zakres:", min(dane$dochod), "-", max(dane$dochod), "\n\n")

Zakres: 40 - 60 

# Statystyki dla frekwencji
cat("FREKWENCJA (%):\n")

FREKWENCJA (%):

cat("Średnia:", mean(dane$frekwencja), "\n")

Średnia: 62.2 

cat("Mediana:", median(dane$frekwencja), "\n")

Mediana: 60 

cat("Odchylenie standardowe:", round(sd(dane$frekwencja), 2), "\n")

Odchylenie standardowe: 10.21 

cat("Zakres:", min(dane$frekwencja), "-", max(dane$frekwencja))

Zakres: 50 - 75

Część 2: Analiza korelacji

# Test korelacji Pearsona
korelacja <- cor.test(dane$dochod, dane$frekwencja)

cat("Współczynnik korelacji (r):", round(korelacja$estimate, 3), "\n")

Współczynnik korelacji (r): 0.994 

cat("P-value:", round(korelacja$p.value, 3), "\n")

P-value: 0.001 

# Interpretacja
if (korelacja$p.value < 0.05) {
  cat("Wynik jest statystycznie istotny (p < 0.05)")
} else {
  cat("Wynik nie jest statystycznie istotny (p ≥ 0.05)")
}

Wynik jest statystycznie istotny (p < 0.05)

Część 3: Model regresji liniowej

# Dopasowanie modelu
model <- lm(frekwencja ~ dochod, data = dane)

# Podstawowe informacje o modelu
summary(model)

Call:
lm(formula = frekwencja ~ dochod, data = dane)

Residuals:
      1       2       3       4       5 
-1.9486  0.3359  0.5100  0.6204  0.4824 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.89647    3.96731  -0.226 0.835748    
dochod       1.25690    0.07822  16.068 0.000524 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.263 on 3 degrees of freedom
Multiple R-squared:  0.9885,    Adjusted R-squared:  0.9847 
F-statistic: 258.2 on 1 and 3 DF,  p-value: 0.0005243

# R-squared - jak dobrze model wyjaśnia dane
cat("\nModel wyjaśnia", round(summary(model)$r.squared * 100, 1), "% zmienności frekwencji")

Model wyjaśnia 98.9 % zmienności frekwencji

Część 4: Wizualizacja

# Wykres rozrzutu z linią regresji
ggplot(dane, aes(x = dochod, y = frekwencja)) +
  geom_point(size = 4, color = "steelblue") +
  geom_smooth(method = "lm", se = TRUE, color = "red", alpha = 0.3) +
  geom_text(aes(label = dzielnica), vjust = -1, size = 3) +
  labs(
    title = "Dochód vs Frekwencja wyborcza",
    subtitle = paste("r =", round(korelacja$estimate, 3)),
    x = "Średni dochód (tys. €)",
    y = "Frekwencja wyborcza (%)"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold"))

`geom_smooth()` using formula = 'y ~ x'

Część 5: Interpretacja wyników

# Pobranie współczynników
wspolczynniki <- coef(model)
przeciecie <- round(wspolczynniki[1], 1)
nachylenie <- round(wspolczynniki[2], 2)

cat("RÓWNANIE REGRESJI:\n")

RÓWNANIE REGRESJI:

cat("Frekwencja =", przeciecie, "+", nachylenie, "× Dochód\n\n")

Frekwencja = -0.9 + 1.26 × Dochód

cat("INTERPRETACJA:\n")

INTERPRETACJA:

cat("• Wzrost dochodu o 1000€ zwiększa frekwencję o", nachylenie, "punktów procentowych\n")

• Wzrost dochodu o 1000€ zwiększa frekwencję o 1.26 punktów procentowych

cat("• Korelacja jest", ifelse(korelacja$estimate > 0, "dodatnia", "ujemna"), 
    "i", ifelse(abs(korelacja$estimate) > 0.7, "silna", 
               ifelse(abs(korelacja$estimate) > 0.5, "umiarkowana", "słaba")))

• Korelacja jest dodatnia i silna

Wnioski

Główne ustalenia:

• Istnieje silna dodatnia korelacja (r = 0.994) między dochodem a frekwencją wyborczą
• Model wyjaśnia 98.9% zmienności w danych
• Dzielnice z wyższymi dochodami mają wyższą frekwencję wyborczą

Ważne ograniczenie:

⚠️ Mała próba (n=5) oznacza, że wyniki należy traktować ostrożnie i nie można ich generalizować na całą populację bez dodatkowych badań.

Praktyczne zastosowanie:

Wyniki sugerują, że działania mające na celu zwiększenie frekwencji wyborczej powinny szczególnie koncentrować się na dzielnicach o niższych dochodach.

Ograniczenia i zastrzeżenia:

⚠️ Krytyczne ograniczenia:

• Bardzo mała próba (n=5) znacznie ogranicza możliwość generalizacji
• Niska moc statystyczna - ryzyko błędów II rodzaju
• Brak kontroli zmiennych zakłócających (wiek, wykształcenie, gęstość zaludnienia)
• Możliwa korelacja pozorna - potrzebne dodatkowe zmienne kontrolne

Rekomendacje dla przyszłych badań:

• Zwiększenie próby do wszystkich dzielnic miasta
• Włączenie zmiennych demograficznych i socjoekonomicznych
• Analiza danych longitudinalnych z kilku cykli wyborczych

10.18 Przykład 2. Związek Między Wielkością Okręgu a Dysproporcjonalnością Wyborczą (1)

Ta analiza bada związek między wielkością okręgu wyborczego (DM) a wskaźnikiem dysproporcjonalności Loosemore-Hanby (LH) w wyborach parlamentarnych w 6 krajach. Indeks Loosemore-Hanby mierzy dysproporcjonalność wyborczą, gdzie wyższe wartości wskazują na większą dysproporcjonalność między głosami a mandatami.

Dane

Warning: package 'knitr' was built under R version 4.4.3

Wielkość Okręgu i Indeks LH według Kraju

Country	DM	LH
A	3	15.50
B	4	14.25
C	5	13.50
D	6	13.50
E	7	13.00
F	8	12.75

Krok 1: Obliczenie wariancji i odchylenia standardowego dla DM i LH

Obliczenia dla Wielkości Okręgu (DM)

Najpierw obliczam średnią wartości DM:

\bar{x}_{DM} = \frac{1}{n}\sum_{i=1}^{n}x_i = \frac{3 + 4 + 5 + 6 + 7 + 8}{6} = \frac{33}{6} = 5.5

Następnie obliczam wariancję używając formuły z korektą Bessela:

\sigma^2_{DM} = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2

Kraj	DM (x_i)	x_i - \bar{x}	(x_i - \bar{x})^2
A	3	3 - 5.5 = -2.5	(-2.5)^2 = 6.25
B	4	4 - 5.5 = -1.5	(-1.5)^2 = 2.25
C	5	5 - 5.5 = -0.5	(-0.5)^2 = 0.25
D	6	6 - 5.5 = 0.5	(0.5)^2 = 0.25
E	7	7 - 5.5 = 1.5	(1.5)^2 = 2.25
F	8	8 - 5.5 = 2.5	(2.5)^2 = 6.25
Suma			17.5

\sigma^2_{DM} = \frac{17.5}{6-1} = \frac{17.5}{5} = 3.5

Odchylenie standardowe to pierwiastek kwadratowy z wariancji:

\sigma_{DM} = \sqrt{\sigma^2_{DM}} = \sqrt{3.5} = 1.871

Obliczenia dla Indeksu LH

Najpierw obliczam średnią wartości LH:

\bar{x}_{LH} = \frac{15.5 + 14.25 + 13.5 + 13.5 + 13 + 12.75}{6} = \frac{82.5}{6} = 13.75

Następnie obliczam wariancję:

Kraj	LH (y_i)	y_i - \bar{y}	(y_i - \bar{y})^2
A	15.5	15.5 - 13.75 = 1.75	(1.75)^2 = 3.0625
B	14.25	14.25 - 13.75 = 0.5	(0.5)^2 = 0.25
C	13.5	13.5 - 13.75 = -0.25	(-0.25)^2 = 0.0625
D	13.5	13.5 - 13.75 = -0.25	(-0.25)^2 = 0.0625
E	13	13 - 13.75 = -0.75	(-0.75)^2 = 0.5625
F	12.75	12.75 - 13.75 = -1	(-1)^2 = 1
Suma			5

\sigma^2_{LH} = \frac{5}{6-1} = \frac{5}{5} = 1

Odchylenie standardowe wynosi:

\sigma_{LH} = \sqrt{\sigma^2_{LH}} = \sqrt{1} = 1

Podsumowanie Zadania 1:

• Wariancja DM (z korektą Bessela): 3.5
• Odchylenie Standardowe DM: 1.871
• Wariancja LH (z korektą Bessela): 1
• Odchylenie Standardowe LH: 1

Krok 2: Obliczenie kowariancji między DM i LH dla tej próby danych

Kowariancja jest obliczana przy użyciu formuły z korektą Bessela:

Cov(DM, LH) = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})

Kraj	DM (x_i)	LH (y_i)	x_i - \bar{x}	y_i - \bar{y}	(x_i - \bar{x})(y_i - \bar{y})
A	3	15.5	-2.5	1.75	(-2.5)(1.75) = -4.375
B	4	14.25	-1.5	0.5	(-1.5)(0.5) = -0.75
C	5	13.5	-0.5	-0.25	(-0.5)(-0.25) = 0.125
D	6	13.5	0.5	-0.25	(0.5)(-0.25) = -0.125
E	7	13	1.5	-0.75	(1.5)(-0.75) = -1.125
F	8	12.75	2.5	-1	(2.5)(-1) = -2.5
Suma					-8.75

Cov(DM, LH) = \frac{-8.75}{5} = -1.75

Kowariancja między DM i LH: -1.75

Ujemna kowariancja wskazuje na odwrotną zależność: gdy wielkość okręgu wzrasta, indeks dysproporcjonalności LH ma tendencję do spadku.

Krok 3: Obliczenie współczynnika korelacji liniowej Pearsona między DM i LH

Współczynnik korelacji Pearsona obliczany jest przy użyciu formuły:

r = \frac{Cov(DM, LH)}{\sigma_{DM} \cdot \sigma_{LH}}

Mamy już obliczone:

• Cov(DM, LH) = -1.75
• \sigma_{DM} = 1.871
• \sigma_{LH} = 1

r = \frac{-1.75}{1.871 \cdot 1} = \frac{-1.75}{1.871} = -0.935

Współczynnik korelacji Pearsona: -0.935

Interpretacja:

Współczynnik korelacji -0.935 wskazuje:

1. Kierunek: Znak ujemny pokazuje odwrotną zależność między wielkością okręgu a indeksem LH.

2. Siła: Wartość bezwzględna 0.935 wskazuje na bardzo silną korelację (blisko -1).

3. Interpretacja praktyczna: Ponieważ wyższe wartości indeksu LH wskazują na większą dysproporcjonalność, ta silna ujemna korelacja sugeruje, że gdy wielkość okręgu wzrasta, dysproporcjonalność wyborcza ma tendencję do znacznego spadku. Innymi słowy, systemy wyborcze z większymi okręgami (więcej przedstawicieli wybieranych z jednego okręgu) zwykle dają bardziej proporcjonalne wyniki (niższa dysproporcjonalność).

Odkrycie to jest zgodne z teorią nauk politycznych, która sugeruje, że większe okręgi zapewniają więcej możliwości mniejszym partiom, aby uzyskać reprezentację, co prowadzi do wyników wyborczych, które lepiej odzwierciedlają rozkład głosów między partiami.

Krok 4: Skonstruowanie modelu regresji liniowej prostej i obliczenie R-kwadrat

Formuła dla regresji liniowej prostej:

LH = \beta_0 + \beta_1 \cdot DM + \varepsilon

Gdzie:

• \beta_0 to wyraz wolny (przecięcie)
• \beta_1 to współczynnik nachylenia dla DM

Formuła do obliczenia \beta_1 to:

\beta_1 = \frac{Cov(DM, LH)}{\sigma^2_{DM}} = \frac{-1.75}{3.5} = -0.5

Aby obliczyć \beta_0, używam:

\beta_0 = \bar{y} - \beta_1 \cdot \bar{x} = 13.75 - (-0.5) \cdot 5.5 = 13.75 + 2.75 = 16.5

Zatem równanie regresji to:

LH = 16.5 - 0.5 \cdot DM

Obliczanie wartości przewidywanych i błędów

Używając naszego równania regresji, obliczam przewidywane wartości LH:

\hat{LH} = 16.5 - 0.5 \cdot DM

Kraj	DM (x_i)	LH (y_i)	Przewidywane LH (\hat{y}_i = 16.5 - 0.5x_i)	Błąd (y_i - \hat{y}_i)	Błąd bezwzględny (|y_i - \hat{y}_i|)	Błąd kwadratowy ([y_i - \hat{y}_i]^2)
A	3	15.5	16.5 - 0.5(3) = 16.5 - 1.5 = 15	15.5 - 15 = 0.5	|0.5| = 0.5	(0.5)^2 = 0.25
B	4	14.25	16.5 - 0.5(4) = 16.5 - 2 = 14.5	14.25 - 14.5 = -0.25	|-0.25| = 0.25	(-0.25)^2 = 0.0625
C	5	13.5	16.5 - 0.5(5) = 16.5 - 2.5 = 14	13.5 - 14 = -0.5	|-0.5| = 0.5	(-0.5)^2 = 0.25
D	6	13.5	16.5 - 0.5(6) = 16.5 - 3 = 13.5	13.5 - 13.5 = 0	|0| = 0	(0)^2 = 0
E	7	13	16.5 - 0.5(7) = 16.5 - 3.5 = 13	13 - 13 = 0	|0| = 0	(0)^2 = 0
F	8	12.75	16.5 - 0.5(8) = 16.5 - 4 = 12.5	12.75 - 12.5 = 0.25	|0.25| = 0.25	(0.25)^2 = 0.0625
Suma					1.5	0.625

Obliczanie R-kwadrat

1. SST (Całkowita suma kwadratów)

SST = \sum_{i=1}^{n}(y_i - \bar{y})^2

Obliczyliśmy już te odchylenia w Kroku 1:

Kraj	LH (y_i)	y_i - \bar{y}	(y_i - \bar{y})^2
A	15.5	1.75	3.0625
B	14.25	0.5	0.25
C	13.5	-0.25	0.0625
D	13.5	-0.25	0.0625
E	13	-0.75	0.5625
F	12.75	-1	1
Suma			5

SST = 5

2. SSR (Suma kwadratów regresji)

SSR = \sum_{i=1}^{n}(\hat{y}_i - \bar{y})^2

Kraj	Przewidywane LH (\hat{y}_i)	\hat{y}_i - \bar{y}	(\hat{y}_i - \bar{y})^2
A	15	15 - 13.75 = 1.25	(1.25)^2 = 1.5625
B	14.5	14.5 - 13.75 = 0.75	(0.75)^2 = 0.5625
C	14	14 - 13.75 = 0.25	(0.25)^2 = 0.0625
D	13.5	13.5 - 13.75 = -0.25	(-0.25)^2 = 0.0625
E	13	13 - 13.75 = -0.75	(-0.75)^2 = 0.5625
F	12.5	12.5 - 13.75 = -1.25	(-1.25)^2 = 1.5625
Suma			4.375

SSR = 4.375

3. SSE (Suma kwadratów błędów)

SSE = \sum_{i=1}^{n}(y_i - \hat{y}_i)^2

Z tabeli powyżej, suma kwadratów błędów wynosi:

SSE = 0.625

4. Weryfikacja

Możemy zweryfikować nasze obliczenia sprawdzając, czy SST = SSR + SSE:

5 = 4.375 + 0.625 = 5 \checkmark

5. Obliczanie R-kwadrat

R^2 = \frac{SSR}{SST} = \frac{4.375}{5} = 0.875

Obliczanie RMSE (Root Mean Square Error)

RMSE jest obliczane przy użyciu formuły:

RMSE = \sqrt{\frac{SSE}{n}} = \sqrt{\frac{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}{n}}

Używając naszego obliczonego SSE:

RMSE = \sqrt{\frac{0.625}{6}} = \sqrt{0.104} \approx 0.323

Obliczanie MAE (Mean Absolute Error)

MAE jest obliczane przy użyciu formuły:

MAE = \frac{\sum_{i=1}^{n}|y_i - \hat{y}_i|}{n}

Używając sum z tabeli:

MAE = \frac{1.5}{6} = 0.25

Model regresji: LH = 16.5 - 0.5 \cdot DM

R-kwadrat: 0.875

RMSE: 0.323

MAE: 0.25

Interpretacja:

1. Równanie regresji: Dla każdego wzrostu jednostkowego wielkości okręgu, indeks dysproporcjonalności LH jest oczekiwany spadek o 0.5 jednostki. Wyraz wolny (16.5) reprezentuje oczekiwany indeks LH, gdy wielkość okręgu wynosi zero (choć nie ma to praktycznego znaczenia, ponieważ wielkość okręgu nie może wynosić zero).

2. R-kwadrat: 0.875 wskazuje, że około 87.5% wariancji w dysproporcjonalności wyborczej (indeks LH) może być wyjaśnione przez wielkość okręgu. Jest to wysoka wartość, sugerująca, że wielkość okręgu jest rzeczywiście silnym predyktorem dysproporcjonalności wyborczej.

3. RMSE i MAE: Niskie wartości RMSE (0.323) i MAE (0.25) wskazują, że model dobrze dopasowuje się do danych, z małymi błędami predykcji.

4. Implikacje polityczne: Odkrycia sugerują, że zwiększanie wielkości okręgu mogłoby być skuteczną strategią reformy wyborczej dla krajów dążących do zmniejszenia dysproporcjonalności między głosami a mandatami. Jednak korzyści marginalne wydają się zmniejszać wraz ze wzrostem wielkości okręgu, jak widać w wzorcu danych.

10.19 Przykład 3. Analiza związku między wielkością okręgu a wskaźnikiem dysproporcjonalności wyborczej (2)

Ta analiza bada związek między wielkością okręgu (DM) a wskaźnikiem dysproporcjonalności Loosemore-Hanby (LH) w wyborach parlamentarnych w 6 krajach. Wskaźnik Loosemore-Hanby mierzy dysproporcjonalność wyborczą, przy czym wyższe wartości wskazują na większą dysproporcjonalność między głosami a mandatami.

Dane

Wielkość okręgu i wskaźnik LH według kraju

Kraj	DM	LH
A	4	12
B	10	8
C	3	15
D	8	10
E	7	6
F	4	13

Krok 1: Obliczenie wariancji i odchylenia standardowego dla DM i LH

Obliczenia dla wielkości okręgu (DM)

Najpierw obliczę średnią wartości DM:

\bar{x}_{DM} = \frac{1}{n}\sum_{i=1}^{n}x_i = \frac{4 + 10 + 3 + 8 + 7 + 4}{6} = \frac{36}{6} = 6

Następnie obliczę wariancję z korektą Bessela, korzystając z wzoru:

s^2_{DM} = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2

Kraj	DM (x_i)	x_i - \bar{x}	(x_i - \bar{x})^2
A	4	4 - 6 = -2	(-2)^2 = 4
B	10	10 - 6 = 4	(4)^2 = 16
C	3	3 - 6 = -3	(-3)^2 = 9
D	8	8 - 6 = 2	(2)^2 = 4
E	7	7 - 6 = 1	(1)^2 = 1
F	4	4 - 6 = -2	(-2)^2 = 4
Suma			38

s^2_{DM} = \frac{38}{6-1} = \frac{38}{5} = 7.6

Odchylenie standardowe to pierwiastek kwadratowy z wariancji:

s_{DM} = \sqrt{s^2_{DM}} = \sqrt{7.6} = 2.757

Obliczenia dla wskaźnika LH

Najpierw obliczę średnią wartości LH:

\bar{y}_{LH} = \frac{12 + 8 + 15 + 10 + 6 + 13}{6} = \frac{64}{6} = 10.667

Następnie obliczę wariancję z korektą Bessela:

Kraj	LH (y_i)	y_i - \bar{y}	(y_i - \bar{y})^2
A	12	12 - 10.667 = 1.333	(1.333)^2 = 1.777
B	8	8 - 10.667 = -2.667	(-2.667)^2 = 7.113
C	15	15 - 10.667 = 4.333	(4.333)^2 = 18.775
D	10	10 - 10.667 = -0.667	(-0.667)^2 = 0.445
E	6	6 - 10.667 = -4.667	(-4.667)^2 = 21.781
F	13	13 - 10.667 = 2.333	(2.333)^2 = 5.443
Suma			55.334

s^2_{LH} = \frac{55.334}{6-1} = \frac{55.334}{5} = 11.067

Odchylenie standardowe to:

s_{LH} = \sqrt{s^2_{LH}} = \sqrt{11.067} = 3.327

Podsumowanie kroku 1:

• Wariancja DM (z korektą Bessela): 7.6
• Odchylenie standardowe DM: 2.757
• Wariancja LH (z korektą Bessela): 11.067
• Odchylenie standardowe LH: 3.327

Krok 2: Obliczenie kowariancji między DM a LH dla tej próby danych

Kowariancja jest obliczana przy użyciu wzoru z korektą Bessela:

Cov(DM, LH) = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})

Kraj	DM (x_i)	LH (y_i)	x_i - \bar{x}	y_i - \bar{y}	(x_i - \bar{x})(y_i - \bar{y})
A	4	12	-2	1.333	(-2)(1.333) = -2.666
B	10	8	4	-2.667	(4)(-2.667) = -10.668
C	3	15	-3	4.333	(-3)(4.333) = -12.999
D	8	10	2	-0.667	(2)(-0.667) = -1.334
E	7	6	1	-4.667	(1)(-4.667) = -4.667
F	4	13	-2	2.333	(-2)(2.333) = -4.666
Suma					-37

Cov(DM, LH) = \frac{-37}{6-1} = \frac{-37}{5} = -7.4

Kowariancja między DM a LH: -7.4

Ujemna kowariancja wskazuje na odwrotną zależność: wraz ze wzrostem wielkości okręgu wskaźnik dysproporcjonalności LH ma tendencję do spadku.

Krok 3: Obliczenie współczynnika korelacji liniowej Pearsona między DM a LH

Współczynnik korelacji Pearsona oblicza się przy użyciu wzoru:

r = \frac{Cov(DM, LH)}{s_{DM} \cdot s_{LH}}

Mamy już obliczone:

• Cov(DM, LH) = -7.4
• s_{DM} = 2.757
• s_{LH} = 3.327

r = \frac{-7.4}{2.757 \cdot 3.327} = \frac{-7.4}{9.172} = -0.807

Współczynnik korelacji Pearsona: -0.807

Interpretacja:

Współczynnik korelacji -0.807 wskazuje:

1. Kierunek: Ujemny znak pokazuje odwrotną zależność między wielkością okręgu a wskaźnikiem LH.

2. Siła: Wartość bezwzględna 0.807 wskazuje na silną korelację (blisko -1).

3. Interpretacja praktyczna: Ponieważ wyższe wartości wskaźnika LH wskazują na większą dysproporcjonalność, ta silna ujemna korelacja sugeruje, że wraz ze wzrostem wielkości okręgu, dysproporcjonalność wyborcza ma tendencję do znacznego spadku. Innymi słowy, systemy wyborcze z większymi okręgami wyborczymi (więcej przedstawicieli wybieranych w okręgu) mają tendencję do generowania bardziej proporcjonalnych wyników (mniejsza dysproporcjonalność).

Ustalenie to jest zgodne z teorią nauk politycznych, która sugeruje, że większe okręgi wyborcze zapewniają mniejszym partiom więcej możliwości uzyskania reprezentacji, co prowadzi do wyników wyborczych, które lepiej odzwierciedlają rozkład głosów między partiami.

Krok 4: Skonstruowanie prostego modelu regresji liniowej i obliczenie R-kwadrat

Użyję wzoru na prostą regresję liniową:

LH = \beta_0 + \beta_1 \cdot DM + \varepsilon

Gdzie:

• \beta_0 to wyraz wolny
• \beta_1 to współczynnik nachylenia dla DM

Wzór na obliczenie \beta_1 z korektą Bessela to:

\beta_1 = \frac{Cov(DM, LH)}{s^2_{DM}} = \frac{-7.4}{7.6} = -0.974

Aby obliczyć \beta_0, użyję:

\beta_0 = \bar{y} - \beta_1 \cdot \bar{x} = 10.667 - (-0.974) \cdot 6 = 10.667 + 5.844 = 16.511

Zatem równanie regresji to:

LH = 16.511 - 0.974 \cdot DM

Obliczanie R-kwadrat

Aby właściwie obliczyć R-kwadrat, muszę obliczyć następujące sumy kwadratów:

• SST (Całkowita suma kwadratów): Mierzy całkowitą zmienność zmiennej zależnej (LH)
• SSR (Suma kwadratów regresji): Mierzy zmienność wyjaśnioną przez model regresji
• SSE (Suma kwadratów błędów): Mierzy niewyjaśnioną zmienność modelu

Najpierw obliczę przewidywane wartości LH, używając naszego równania regresji:

\hat{LH} = 16.511 - 0.974 \cdot DM

Kraj	DM (x_i)	LH (y_i)	Przewidywane LH (\hat{y}_i = 16.511 - 0.974x_i)
A	4	12	16.511 - 0.974(4) = 16.511 - 3.896 = 12.615
B	10	8	16.511 - 0.974(10) = 16.511 - 9.74 = 6.771
C	3	15	16.511 - 0.974(3) = 16.511 - 2.922 = 13.589
D	8	10	16.511 - 0.974(8) = 16.511 - 7.792 = 8.719
E	7	6	16.511 - 0.974(7) = 16.511 - 6.818 = 9.693
F	4	13	16.511 - 0.974(4) = 16.511 - 3.896 = 12.615

Teraz, obliczę każdą sumę kwadratów:

1. SST (Całkowita suma kwadratów)

SST = \sum_{i=1}^{n}(y_i - \bar{y})^2

Już obliczyliśmy te odchylenia w Kroku 1:

Kraj	LH (y_i)	y_i - \bar{y}	(y_i - \bar{y})^2
A	12	1.333	1.777
B	8	-2.667	7.113
C	15	4.333	18.775
D	10	-0.667	0.445
E	6	-4.667	21.781
F	13	2.333	5.443
Suma			55.334

SST = 55.334

2. SSR (Suma kwadratów regresji)

SSR = \sum_{i=1}^{n}(\hat{y}_i - \bar{y})^2

Kraj	Przewidywane LH (\hat{y}_i)	\hat{y}_i - \bar{y}	(\hat{y}_i - \bar{y})^2
A	12.615	12.615 - 10.667 = 1.948	(1.948)^2 = 3.795
B	6.771	6.771 - 10.667 = -3.896	(-3.896)^2 = 15.178
C	13.589	13.589 - 10.667 = 2.922	(2.922)^2 = 8.538
D	8.719	8.719 - 10.667 = -1.948	(-1.948)^2 = 3.795
E	9.693	9.693 - 10.667 = -0.974	(-0.974)^2 = 0.949
F	12.615	12.615 - 10.667 = 1.948	(1.948)^2 = 3.795
Suma			36.05

SSR = 36.05

3. SSE (Suma kwadratów błędów)

SSE = \sum_{i=1}^{n}(y_i - \hat{y}_i)^2

Kraj	LH (y_i)	Przewidywane LH (\hat{y}_i)	y_i - \hat{y}_i	(y_i - \hat{y}_i)^2
A	12	12.615	12 - 12.615 = -0.615	(-0.615)^2 = 0.378
B	8	6.771	8 - 6.771 = 1.229	(1.229)^2 = 1.510
C	15	13.589	15 - 13.589 = 1.411	(1.411)^2 = 1.991
D	10	8.719	10 - 8.719 = 1.281	(1.281)^2 = 1.641
E	6	9.693	6 - 9.693 = -3.693	(-3.693)^2 = 13.638
F	13	12.615	13 - 12.615 = 0.385	(0.385)^2 = 0.148
Suma				19.306

SSE = 19.306

4. Weryfikacja

Możemy zweryfikować nasze obliczenia, sprawdzając czy SST = SSR + SSE:

55.334 \approx 36.05 + 19.306 = 55.356

Niewielka różnica (0.022) wynika z zaokrągleń w obliczeniach.

5. Obliczanie R-kwadrat

R^2 = \frac{SSR}{SST} = \frac{36.05}{55.334} = 0.652

Alternatywnie, możemy też obliczyć:

R^2 = 1 - \frac{SSE}{SST} = 1 - \frac{19.306}{55.334} = 1 - 0.349 = 0.651

Drobna różnica wynika z zaokrągleń.

Obliczanie RMSE (Pierwiastek średniego błędu kwadratowego)

Obliczanie RMSE (Pierwiastek średniego błędu kwadratowego)

RMSE oblicza się przy użyciu wzoru:

RMSE = \sqrt{\frac{SSE}{n}} = \sqrt{\frac{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}{n}}

Korzystając z naszej obliczonej SSE:

RMSE = \sqrt{\frac{19.306}{6}} = \sqrt{3.218} = 1.794

Korekta Bessela (dzielenie przez n-1 zamiast n) stosuje się do estymacji wariancji próby, ale nie jest standardowo stosowana przy obliczaniu RMSE, gdyż RMSE jest miarą błędu predykcji, a nie estymatorem parametru populacji.

Obliczanie MAE (Średni błąd bezwzględny)

MAE oblicza się przy użyciu wzoru:

MAE = \frac{1}{n}\sum_{i=1}^{n}|y_i - \hat{y}_i|

Kraj	LH (y_i)	Przewidywane LH (\hat{y}_i)	|y_i - \hat{y}_i|
A	12	12.615	|12 - 12.615| = 0.615
B	8	6.771	|8 - 6.771| = 1.229
C	15	13.589	|15 - 13.589| = 1.411
D	10	8.719	|10 - 8.719| = 1.281
E	6	9.693	|6 - 9.693| = 3.693
F	13	12.615	|13 - 12.615| = 0.385
Suma			8.614

MAE = \frac{8.614}{6} = 1.436

Model regresji: LH = 16.511 - 0.974 \cdot DM

R-kwadrat: 0.651

RMSE: 1.794

MAE: 1.436

Interpretacja:

1. Równanie regresji: Dla każdego jednostkowego wzrostu wielkości okręgu, wskaźnik dysproporcjonalności LH zmniejsza się o 0.974 jednostki. Wyraz wolny (16.511) reprezentuje oczekiwany wskaźnik LH, gdy wielkość okręgu wynosi zero (choć nie ma to praktycznego znaczenia, ponieważ wielkość okręgu nie może wynosić zero).

2. R-kwadrat: 0.651 wskazuje, że około 65.1% wariancji dysproporcjonalności wyborczej (wskaźnik LH) może być wyjaśnione przez wielkość okręgu. Jest to dość wysoka wartość, sugerująca, że wielkość okręgu jest rzeczywiście silnym predyktorem dysproporcjonalności wyborczej, choć mniejszym niż w poprzednim zestawie danych.

3. RMSE: Wartość 1.794 informuje nas o przeciętnym błędzie prognozy modelu. Jest to miara dokładności przewidywań modelu wyrażona w jednostkach zmiennej zależnej (LH).

4. MAE: Wartość 1.436 informuje nas o przeciętnym bezwzględnym błędzie prognozy modelu. W porównaniu z RMSE, MAE jest mniej czuły na wartości odstające, co potwierdza, że niektóre obserwacje (np. dla kraju E) mają stosunkowo duży błąd predykcji.

5. Implikacje polityczne: Wyniki sugerują, że zwiększenie wielkości okręgu mogłoby być skuteczną strategią reform wyborczych dla krajów starających się zmniejszyć dysproporcjonalność między głosami a mandatami. Jednakże, korzyści marginalne wydają się zmniejszać wraz ze wzrostem wielkości okręgu, jak widać we wzorcu danych.

10.20 Przykład 4. Lęk vs. Wyniki Egzaminu: Analiza Korelacji i Regresji

W tym tutorialu zbadamy związek między poziomem lęku przed egzaminem a wynikami egzaminacyjnymi wśród studentów uniwersytetu. Badania sugerują, że podczas gdy niewielka ilość lęku może być motywująca, nadmierny lęk zazwyczaj pogarsza wyniki poprzez zmniejszoną koncentrację, zakłócenia pamięci roboczej i objawy fizyczne (Yerkes-Dodson law). Przeanalizujemy dane od 8 studentów, aby zrozumieć ten związek matematycznie.

Prezentacja Danych

Zbiór Danych

Zebraliśmy dane od 8 studentów, mierząc:

• X: Wynik lęku przed testem (skala 1-10, gdzie 1 = bardzo niski, 10 = bardzo wysoki)
• Y: Wyniki egzaminu (wynik procentowy)

# Nasze dane
lęk <- c(2.5, 3.2, 4.1, 4.8, 5.6, 6.3, 7.0, 7.9)      # Wyniki lęku
wyniki <- c(80, 85, 78, 82, 77, 74, 68, 72)            # Wyniki egzaminu (%)

# Stworzenie ramki danych dla łatwego przeglądu
dane <- data.frame(
  Student = 1:8,
  Lęk = lęk,
  Wyniki = wyniki
)
print(dane)

  Student Lęk Wyniki
1       1 2.5     80
2       2 3.2     85
3       3 4.1     78
4       4 4.8     82
5       5 5.6     77
6       6 6.3     74
7       7 7.0     68
8       8 7.9     72

Wstępna Wizualizacja

Najpierw zwizualizujmy nasze dane, aby uzyskać intuicyjne zrozumienie związku:

library(ggplot2)

# Stworzenie wykresu punktowego
ggplot(dane, aes(x = Lęk, y = Wyniki)) +
  geom_point(size = 4, color = "darkblue") +
  labs(
    title = "Lęk przed Testem vs. Wyniki Egzaminu",
    x = "Wynik Lęku przed Testem",
    y = "Wyniki Egzaminu (%)"
  ) +
  theme_minimal() +
  scale_x_continuous(breaks = seq(2, 8, 1)) +
  scale_y_continuous(breaks = seq(65, 90, 5))

Statystyki Opisowe

Obliczanie Średnich

Średnia to wartość przeciętna, obliczana przez zsumowanie wszystkich obserwacji i podzielenie przez liczbę obserwacji.

Średnia Wyników Lęku (X): \bar{X} = \frac{\sum_{i=1}^{n} X_i}{n} = \frac{2,5 + 3,2 + 4,1 + 4,8 + 5,6 + 6,3 + 7,0 + 7,9}{8}

\bar{X} = \frac{41,4}{8} = 5,175

Średnia Wyników Egzaminu (Y): \bar{Y} = \frac{\sum_{i=1}^{n} Y_i}{n} = \frac{80 + 85 + 78 + 82 + 77 + 74 + 68 + 72}{8}

\bar{Y} = \frac{616}{8} = 77

# Weryfikacja naszych obliczeń
średnia_x <- mean(lęk)
średnia_y <- mean(wyniki)
cat("Średni Lęk:", średnia_x, "\n")

Średni Lęk: 5.175 

cat("Średnie Wyniki:", średnia_y, "\n")

Średnie Wyniki: 77 

Obliczanie Wariancji i Odchylenia Standardowego

Wariancja mierzy, jak bardzo rozproszone są dane od średniej. Używamy wzoru na wariancję próbkową (dzieląc przez n-1).

Wariancja X:

Najpierw obliczamy odchylenia od średniej (X_i - \bar{X}):

Student	X_i	X_i - \bar{X}	(X_i - \bar{X})^2
1	2,5	2,5 - 5,175 = -2,675	(-2,675)^2 = 7,1556
2	3,2	3,2 - 5,175 = -1,975	(-1,975)^2 = 3,9006
3	4,1	4,1 - 5,175 = -1,075	(-1,075)^2 = 1,1556
4	4,8	4,8 - 5,175 = -0,375	(-0,375)^2 = 0,1406
5	5,6	5,6 - 5,175 = 0,425	(0,425)^2 = 0,1806
6	6,3	6,3 - 5,175 = 1,125	(1,125)^2 = 1,2656
7	7,0	7,0 - 5,175 = 1,825	(1,825)^2 = 3,3306
8	7,9	7,9 - 5,175 = 2,725	(2,725)^2 = 7,4256
Suma			24,555

s^2_X = \frac{\sum(X_i - \bar{X})^2}{n-1} = \frac{24,555}{7} = 3,5079

s_X = \sqrt{3,5079} = 1,8730

Wariancja Y:

Student	Y_i	Y_i - \bar{Y}	(Y_i - \bar{Y})^2
1	80	80 - 77 = 3	(3)^2 = 9
2	85	85 - 77 = 8	(8)^2 = 64
3	78	78 - 77 = 1	(1)^2 = 1
4	82	82 - 77 = 5	(5)^2 = 25
5	77	77 - 77 = 0	(0)^2 = 0
6	74	74 - 77 = -3	(-3)^2 = 9
7	68	68 - 77 = -9	(-9)^2 = 81
8	72	72 - 77 = -5	(-5)^2 = 25
Suma			214

s^2_Y = \frac{\sum(Y_i - \bar{Y})^2}{n-1} = \frac{214}{7} = 30,5714

s_Y = \sqrt{30,5714} = 5,5291

# Weryfikacja wariancji i odchylenia standardowego
cat("Wariancja Lęku:", var(lęk), "\n")

Wariancja Lęku: 3.507857 

cat("Odch. Stand. Lęku:", sd(lęk), "\n")

Odch. Stand. Lęku: 1.872927 

cat("Wariancja Wyników:", var(wyniki), "\n")

Wariancja Wyników: 30.57143 

cat("Odch. Stand. Wyników:", sd(wyniki), "\n")

Odch. Stand. Wyników: 5.529144 

Kowariancja i Korelacja

Obliczanie Kowariancji

Kowariancja mierzy, jak dwie zmienne zmieniają się razem. Ujemna kowariancja wskazuje, że gdy jedna zmienna wzrasta, druga ma tendencję do spadku.

s_{XY} = \frac{\sum_{i=1}^{n}(X_i - \bar{X})(Y_i - \bar{Y})}{n-1}

Obliczmy iloczyny dla każdej obserwacji:

Student	(X_i - \bar{X})	(Y_i - \bar{Y})	(X_i - \bar{X})(Y_i - \bar{Y})
1	-2,675	3	(-2,675) \times (3) = -8,025
2	-1,975	8	(-1,975) \times (8) = -15,800
3	-1,075	1	(-1,075) \times (1) = -1,075
4	-0,375	5	(-0,375) \times (5) = -1,875
5	0,425	0	(0,425) \times (0) = 0
6	1,125	-3	(1,125) \times (-3) = -3,375
7	1,825	-9	(1,825) \times (-9) = -16,425
8	2,725	-5	(2,725) \times (-5) = -13,625
Suma			-60,2

s_{XY} = \frac{-60,2}{7} = -8,6

Obliczanie Współczynnika Korelacji Pearsona

Współczynnik korelacji standaryzuje kowariancję do skali od -1 do +1:

r = \frac{s_{XY}}{s_X \cdot s_Y} = \frac{-8,6}{1,8730 \times 5,5291} = \frac{-8,6}{10,3560} = -0,831

To daje nam korelację -0,831, wskazującą na silny ujemny związek między lękiem a wynikami.

# Weryfikacja korelacji
rzeczywista_kor <- cor(lęk, wyniki)
cat("Współczynnik korelacji Pearsona:", rzeczywista_kor, "\n")

Współczynnik korelacji Pearsona: -0.8304618 

Prosta Regresja MNK

Problem Modelowania Związków

Kiedy obserwujemy związek między dwiema zmiennymi, chcemy znaleźć model matematyczny, który:

1. Opisuje związek
2. Pozwala nam dokonywać prognoz
3. Kwantyfikuje siłę związku

Najprostszym modelem jest linia prosta: Y = \beta_0 + \beta_1 X + \epsilon

Gdzie:

• Y to zmienna wynikowa (wyniki)
• X to zmienna predykcyjna (lęk)
• \beta_0 to punkt przecięcia (wyniki gdy lęk = 0)
• \beta_1 to nachylenie (zmiana wyników na jednostkę zmiany lęku)
• \epsilon to składnik błędu (niewyjaśniona zmienność)

Idea Sumy Kwadratów Błędów (SSE)

Wyobraź sobie próbę narysowania linii przez nasze punkty danych. Jest nieskończenie wiele linii, które moglibyśmy narysować! Niektóre przeszłyby przez środek punktów, inne mogłyby być za wysokie lub za niskie, za strome lub za płaskie. Potrzebujemy systematycznego sposobu określenia, która linia jest “najlepsza”.

Czym są Błędy (Reszty)?

Dla każdej linii, którą narysujemy, każdy punkt danych będzie miał błąd lub resztę - pionową odległość od punktu do linii. To reprezentuje, jak “błędna” jest nasza prognoza dla tego punktu.

• Błąd dodatni: Rzeczywista wartość jest powyżej przewidywanej wartości (niedoszacowaliśmy)
• Błąd ujemny: Rzeczywista wartość jest poniżej przewidywanej wartości (przeszacowaliśmy)

Dlaczego Podnosimy Błędy do Kwadratu?

Proste dodawanie błędów nie zadziała, ponieważ błędy dodatnie i ujemne się znoszą. Moglibyśmy użyć wartości bezwzględnych, ale podnoszenie do kwadratu ma kilka zalet:

1. Wygoda matematyczna: Funkcje kwadratowe są różniczkowalne, co ułatwia znalezienie minimum za pomocą rachunku różniczkowego
2. Penalizuje większe błędy bardziej: Kilka dużych błędów jest gorsze niż wiele małych błędów
3. Tworzy gładką, miskowatą funkcję: To gwarantuje unikalne minimum

Wzór SSE

Suma Kwadratów Błędów to: SSE = \sum_{i=1}^{n}(Y_i - \hat{Y_i})^2 = \sum_{i=1}^{n}(Y_i - (\beta_0 + \beta_1 X_i))^2

Znajdowanie Najlepszej Linii

“Najlepsza” linia to ta, która minimalizuje SSE. Używając rachunku różniczkowego (biorąc pochodne względem \beta_0 i \beta_1 i przyrównując je do zera), otrzymujemy wzory MNK.

Estymatory MNK

Metoda Najmniejszych Kwadratów (MNK) znajduje wartości \beta_0 i \beta_1, które minimalizują SSE:

Estymator nachylenia: \hat{\beta_1} = \frac{s_{XY}}{s^2_X} = \frac{\sum(X_i - \bar{X})(Y_i - \bar{Y})}{\sum(X_i - \bar{X})^2}

Estymator punktu przecięcia: \hat{\beta_0} = \bar{Y} - \hat{\beta_1}\bar{X}

Obliczanie Parametrów MNK

Używając naszych obliczonych wartości:

• s_{XY} = -8,6
• s^2_X = 3,5079
• \bar{X} = 5,175
• \bar{Y} = 77

Krok 1: Oblicz nachylenie \hat{\beta_1} = \frac{s_{XY}}{s^2_X} = \frac{-8,6}{3,5079} = -2,451

Krok 2: Oblicz punkt przecięcia \hat{\beta_0} = \bar{Y} - \hat{\beta_1}\bar{X} = 77 - (-2,451 \times 5,175) = 77 + 12,684 = 89,684

# Weryfikacja z R
model <- lm(wyniki ~ lęk)
coef(model)

(Intercept)         lęk 
  89.687233   -2.451639 

Nasze równanie regresji to: \hat{Y} = 89,684 - 2,451X

To oznacza:

• Gdy lęk = 0, przewidywane wyniki = 89,68%
• Za każdy wzrost lęku o 1 punkt, wyniki spadają o 2,45%

Wykres Modelu

library(ggplot2)

# Stworzenie kompleksowego wykresu
ggplot(data.frame(lęk, wyniki), aes(x = lęk, y = wyniki)) +
  geom_point(size = 4, color = "darkblue") +
  geom_smooth(method = "lm", se = FALSE, color = "red", linewidth = 1) +
  geom_text(aes(label = paste0("(", round(lęk, 1), ", ", wyniki, ")")),
            vjust = -1, size = 3) +
  annotate("text", x = 3, y = 70, 
           label = paste0("ŷ = ", round(coef(model)[1], 2), " - ", 
                         abs(round(coef(model)[2], 2)), "x"),
           size = 5, color = "red") +
  labs(
    title = "Linia Regresji: Wyniki vs. Lęk",
    subtitle = "Wyższy lęk jest związany z niższymi wynikami egzaminu",
    x = "Wynik Lęku przed Testem",
    y = "Wyniki Egzaminu (%)"
  ) +
  theme_minimal() +
  scale_x_continuous(breaks = seq(2, 8, 1)) +
  scale_y_continuous(breaks = seq(65, 90, 5))

`geom_smooth()` using formula = 'y ~ x'

Ocena Modelu

Dekompozycja Wariancji

Całkowita zmienność w Y może być rozłożona na dwie części:

SST = SSE + SSR

Gdzie:

• SST (Sum of Squares Total): Całkowita zmienność w Y
• SSE (Sum of Squared Errors): Niewyjaśniona zmienność
  
• SSR (Sum of Squares Regression): Wyjaśniona zmienność

Tip

Wyobraź sobie, że chcesz zrozumieć, dlaczego pensje w firmie się różnią. Jedni zarabiają 40 000, inni 80 000, a jeszcze inni 120 000. Mamy więc zmienność wynagrodzeń — nie są takie same.

Całkowita zmienność (SST)

To pytanie: „Jak bardzo wszystkie pensje są rozproszone wokół średniej pensji?” Jeśli średnia to 70 000, to SST mierzy, o ile każda pensja różni się od 70 000, podnosi te różnice do kwadratu (żeby były dodatnie) i sumuje. To całkowita ilość zmienności, którą próbujemy wyjaśnić.

Zmienność wyjaśniona (SSR)

Załóżmy, że budujemy model przewidujący pensję na podstawie lat doświadczenia. Model może mówić:

• 2 lata doświadczenia → przewiduje 50 000
• 5 lat doświadczenia → przewiduje 70 000
• 10 lat doświadczenia → przewiduje 100 000

SSR mierzy, jak bardzo te przewidywania różnią się od średniej. To ta część zmienności, którą model „wyjaśnia” relacją z doświadczeniem. Innymi słowy: „tę część różnic w pensjach da się przypisać różnym poziomom doświadczenia”.

Zmienność niewyjaśniona (SSE)

To to, co zostaje — część, której model nie tłumaczy. Może być tak, że dwie osoby mają po 5 lat doświadczenia, ale jedna zarabia 65 000, a druga 75 000. Model obu przewidział 70 000. Te różnice względem przewidywań (błędy) reprezentują zmienność wynikającą z innych czynników — np. edukacja, wyniki pracy, umiejętności negocjacyjne albo po prostu losowość.

10.21 Kluczowa intuicja

Piękne jest to, że te trzy wielkości zawsze spełniają zależność: Całkowita zmienność = Zmienność wyjaśniona + Zmienność niewyjaśniona czyli SST = SSR + SSE.

Możesz myśleć o tym jak o „wykresie kołowym powodów, dlaczego pensje się różnią”:

• Jeden kawałek to „różnice wyjaśnione doświadczeniem” (SSR),
• Drugi kawałek to „różnice z innych powodów” (SSE),
• Razem dają całość (SST).

10.22 Dlaczego to ważne

To rozbicie pozwala policzyć współczynnik determinacji R^2:

R^2 = \frac{\mathrm{SSR}}{\mathrm{SST}} = \frac{\text{Zmienność wyjaśniona}}{\text{Całkowita zmienność}}.

Jeśli R^2 = 0{,}70, to model wyjaśnia 70% tego, dlaczego wartości Y (tu: pensje) różnią się między sobą. Pozostałe 30% to czynniki nieuwzględnione lub szum losowy.

Pomyśl o tym jak o rozwiązywaniu zagadki: SST to cała zagadka do rozwiązania, SSR to część już rozwiązana, a SSE to to, co wciąż pozostaje do wyjaśnienia!

Obliczmy każdą:

Krok 1: Oblicz przewidywane wartości

Używając \hat{Y} = 89,684 - 2,451X:

i	X_i	Y_i	\hat{Y}_i = 89,684 - 2,451X_i
1	2,5	80	89,684 - 2,451(2,5) = 83,556
2	3,2	85	89,684 - 2,451(3,2) = 81,841
3	4,1	78	89,684 - 2,451(4,1) = 79,635
4	4,8	82	89,684 - 2,451(4,8) = 77,919
5	5,6	77	89,684 - 2,451(5,6) = 75,968
6	6,3	74	89,684 - 2,451(6,3) = 74,253
7	7,0	68	89,684 - 2,451(7,0) = 72,527
8	7,9	72	89,684 - 2,451(7,9) = 70,321

Krok 2: Oblicz sumy kwadratów

SST (Całkowita zmienność): SST = \sum(Y_i - \bar{Y})^2

Z naszych wcześniejszych obliczeń wariancji: SST = (n-1) \times s^2_Y = 7 \times 30,5714 = 214

SSE (Niewyjaśniona zmienność): SSE = \sum(Y_i - \hat{Y}_i)^2

i	Y_i	\hat{Y}_i	Y_i - \hat{Y}_i	(Y_i - \hat{Y}_i)^2
1	80	83,556	-3,556	12,645
2	85	81,841	3,159	9,979
3	78	79,635	-1,635	2,673
4	82	77,919	4,081	16,655
5	77	75,968	1,032	1,065
6	74	74,253	-0,253	0,064
7	68	72,527	-4,527	20,494
8	72	70,321	1,679	2,819
Suma				66,394

SSR (Wyjaśniona zmienność): SSR = SST - SSE = 214 - 66,394 = 147,606

# Weryfikacja obliczeń
tabela_anova <- anova(model)
cat("SSR (Regresja):", tabela_anova$`Sum Sq`[1], "\n")

SSR (Regresja): 147.5887 

cat("SSE (Reszty):", tabela_anova$`Sum Sq`[2], "\n")

SSE (Reszty): 66.41132 

cat("SST (Całkowita):", sum(tabela_anova$`Sum Sq`), "\n")

SST (Całkowita): 214 

R-kwadrat (Współczynnik Determinacji)

R-kwadrat mówi nam, jaka część całkowitej zmienności w Y jest wyjaśniona przez nasz model:

R^2 = \frac{SSR}{SST} = \frac{147,606}{214} = 0,690

To oznacza, że nasz model wyjaśnia 69,0% zmienności w wynikach egzaminu.

Alternatywny wzór używający korelacji: R^2 = r^2 = (-0,831)^2 = 0,691

# Weryfikacja R-kwadrat
cat("R-kwadrat:", summary(model)$r.squared, "\n")

R-kwadrat: 0.6896667 

cat("Korelacja do kwadratu:", cor(lęk, wyniki)^2, "\n")

Korelacja do kwadratu: 0.6896667 

Wielkość Efektu

Współczynnik nachylenia \hat{\beta_1} = -2,451 to nasza wielkość efektu. Mówi nam:

• Wielkość: Każdy wzrost lęku o 1 punkt jest związany z 2,45% spadkiem wyników
• Znaczenie praktyczne: Student przechodzący od niskiego lęku (3) do wysokiego lęku (7) doświadczyłby oczekiwanego spadku wyników o 2,451 \times 4 = 9,80\%

Standaryzowana wielkość efektu (współczynnik korelacji): Korelacja r = -0,831 wskazuje na silny ujemny związek.

Wytyczne Cohena dla interpretacji korelacji:

• Mały efekt: |r| = 0,10
• Średni efekt: |r| = 0,30
• Duży efekt: |r| = 0,50

Nasze |r| = 0,831 reprezentuje dużą wielkość efektu.

Przedziały Ufności i Istotność Statystyczna

Błąd Standardowy Regresji

Najpierw potrzebujemy błędu standardowego reszt:

s_e = \sqrt{\frac{SSE}{n-2}} = \sqrt{\frac{66,394}{6}} = \sqrt{11,066} = 3,327

Błąd Standardowy Nachylenia

Błąd standardowy \hat{\beta_1} to:

BE(\hat{\beta_1}) = \frac{s_e}{\sqrt{\sum(X_i - \bar{X})^2}} = \frac{3,327}{\sqrt{24,555}} = \frac{3,327}{4,955} = 0,671

95% Przedział Ufności dla Nachylenia

Mówiąc prościej: Przedział ufności daje nam zakres prawdopodobnych wartości dla naszego prawdziwego nachylenia. Gdybyśmy powtórzyli to badanie wiele razy, 95% obliczonych przez nas przedziałów zawierałoby prawdziwą wartość nachylenia.

Wzór używa wartości krytycznej (około 2,45 dla 6 stopni swobody):

PU = \hat{\beta_1} \pm (wartość\_krytyczna \times BE(\hat{\beta_1})) PU = -2,451 \pm (2,45 \times 0,671) PU = -2,451 \pm 1,644 PU = [-4,095, -0,807]

Interpretacja: Jesteśmy w 95% pewni, że prawdziwa zmiana wyników na jednostkę zmiany lęku mieści się między -4,10% a -0,81%.

Istotność Statystyczna

Aby sprawdzić, czy związek jest statystycznie istotny (tj. nie jest przypadkowy), obliczamy statystykę t:

t = \frac{\hat{\beta_1}}{BE(\hat{\beta_1})} = \frac{-2,451}{0,671} = -3,653

Mówiąc prościej: Ta wartość t mówi nam, o ile błędów standardowych nasze nachylenie odbiega od zera. Wartość bezwzględna 3,65 jest dość duża (zazwyczaj wartości przekraczające ±2,45 są uważane za istotne dla naszej wielkości próby), dostarczając silnego dowodu na rzeczywisty ujemny związek między lękiem a wynikami.

# Weryfikacja obliczeń z R
summary(model)

Call:
lm(formula = wyniki ~ lęk)

Residuals:
   Min     1Q Median     3Q    Max 
-4.526 -2.116  0.400  2.050  4.081 

Coefficients:
            Estimate Std. Error t value    Pr(>|t|)    
(Intercept)  89.6872     3.6682  24.450 0.000000308 ***
lęk          -2.4516     0.6714  -3.652      0.0107 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.327 on 6 degrees of freedom
Multiple R-squared:  0.6897,    Adjusted R-squared:  0.6379 
F-statistic: 13.33 on 1 and 6 DF,  p-value: 0.01069

confint(model)

                2.5 %     97.5 %
(Intercept) 80.711582 98.6628838
lęk         -4.094474 -0.8088048

Podsumowanie i Interpretacja

Czego się Nauczyliśmy

1. Istnieje ujemny związek między lękiem przed testem a wynikami egzaminu (r = -0,831)
2. Związek jest umiarkowanie silny - lęk wyjaśnia 69,0% zmienności w wynikach
3. Efekt jest znaczący - każdy wzrost lęku o 1 punkt jest związany z około 2,5% spadkiem wyników
4. Związek jest statystycznie istotny - bardzo mało prawdopodobne, żeby był przypadkowy

Implikacje Praktyczne

Nasza analiza sugeruje, że wysoki lęk przed testem pogarsza wyniki, prawdopodobnie poprzez:

• Zakłócenia poznawcze (niepokojące myśli konkurują o pamięć roboczą)
• Objawy fizyczne (pocenie się, szybkie bicie serca), które odwracają uwagę od zadania
• Negatywny monolog wewnętrzny zmniejszający pewność siebie i motywację
• Zachowania podczas testu (pośpiech, podważanie odpowiedzi)

Rekomendacje na Podstawie Ustaleń

Biorąc pod uwagę silny ujemny związek, interwencje mogłyby obejmować:

• Nauczanie technik zarządzania lękiem (głębokie oddychanie, progresywne rozluźnianie mięśni)
• Restrukturyzację poznawczą w celu radzenia sobie z katastroficznym myśleniem
• Trening umiejętności nauki w celu zwiększenia pewności siebie w przygotowaniu
• Testy próbne w celu zmniejszenia strachu przed nieznanym

Ograniczenia Naszej Analizy

1. Mała wielkość próby (n = 8) ogranicza możliwość uogólnienia
2. Założenie liniowości - związek może być krzywoliniowy (może istnieć optymalny lęk)
3. Inne zmienne nieuwzględnione (czas przygotowania, zdolności, sen, itp.)
4. Przyczynowość vs. korelacja - nie możemy udowodnić, że lęk powoduje słabe wyniki
5. Samoopisowy lęk - subiektywne miary mogą nie odzwierciedlać pobudzenia fizjologicznego

Dodatek: Kompletny Kod R

# Kompletny kod analizy
lęk <- c(2.5, 3.2, 4.1, 4.8, 5.6, 6.3, 7.0, 7.9)
wyniki <- c(80, 85, 78, 82, 77, 74, 68, 72)

# Statystyki opisowe
mean(lęk); sd(lęk)

[1] 5.175

[1] 1.872927

mean(wyniki); sd(wyniki)

[1] 77

[1] 5.529144

# Korelacja
cor(lęk, wyniki)

[1] -0.8304618

# Regresja
model <- lm(wyniki ~ lęk)
summary(model)

Call:
lm(formula = wyniki ~ lęk)

Residuals:
   Min     1Q Median     3Q    Max 
-4.526 -2.116  0.400  2.050  4.081 

Coefficients:
            Estimate Std. Error t value    Pr(>|t|)    
(Intercept)  89.6872     3.6682  24.450 0.000000308 ***
lęk          -2.4516     0.6714  -3.652      0.0107 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.327 on 6 degrees of freedom
Multiple R-squared:  0.6897,    Adjusted R-squared:  0.6379 
F-statistic: 13.33 on 1 and 6 DF,  p-value: 0.01069

confint(model)

                2.5 %     97.5 %
(Intercept) 80.711582 98.6628838
lęk         -4.094474 -0.8088048

# Wizualizacja
library(ggplot2)
ggplot(data.frame(lęk, wyniki), aes(x = lęk, y = wyniki)) +
  geom_point(size = 4) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

# Diagnostyka
plot(model)

# Tabela ANOVA
anova(model)

Analysis of Variance Table

Response: wyniki
          Df  Sum Sq Mean Sq F value  Pr(>F)  
lęk        1 147.589 147.589  13.334 0.01069 *
Residuals  6  66.411  11.069                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

---

10.23 Przykład 5. District Magnitude and Electoral Disproportionality: A Comparative Analysis (*) [EN]

Data Generating Process

Let’s set up a DGP where:

\begin{aligned} & Y_{\text{Gallagher}} = 12 - 0.8X_{\text{DM}} + \varepsilon \\ & \varepsilon \sim \mathcal{N}(0, 1) \\ & X_{\text{DM}} \in \{3, 5, 7, 10, 12, 15\} \end{aligned}

# DGP
magnitude <- c(3, 5, 7, 10, 12, 15)
epsilon <- rnorm(6, mean = 0, sd = 1)
gallagher <- 12 - 0.8 * magnitude + epsilon

# Round (sampled from the DGP) Gallagher indices to one decimal place
gallagher <- round(c(9.0, 7.8, 9.2, 4.1, 2.5, 1.7), 1)

# Show data
data.frame(
  District_Magnitude = magnitude,
  Gallagher_Index = gallagher
)

  District_Magnitude Gallagher_Index
1                  3             9.0
2                  5             7.8
3                  7             9.2
4                 10             4.1
5                 12             2.5
6                 15             1.7

# Initial model check
model <- lm(gallagher ~ magnitude)
summary(model)

Call:
lm(formula = gallagher ~ magnitude)

Residuals:
      1       2       3       4       5       6 
-0.6516 -0.4628  2.3260 -0.6908 -0.9020  0.3813 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  11.7349     1.3034   9.003 0.000843 ***
magnitude    -0.6944     0.1359  -5.110 0.006934 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.368 on 4 degrees of freedom
Multiple R-squared:  0.8672,    Adjusted R-squared:  0.834 
F-statistic: 26.11 on 1 and 4 DF,  p-value: 0.006934

Descriptive Statistics

Means: \bar{X} = \frac{3 + 5 + 7 + 10 + 12 + 15}{6} = \frac{52}{6} = 8.6667

\bar{Y} = \frac{9.0 + 7.8 + 9.2 + 4.1 + 2.5 + 1.7}{6} = \frac{34.3}{6} = 5.7167

# Verification
mean(magnitude)

[1] 8.666667

mean(gallagher)

[1] 5.716667

Variances: s^2_X = \frac{\sum(X_i - \bar{X})^2}{n-1}

Deviations for X:

• (3 - 8.6667) = -5.6667
• (5 - 8.6667) = -3.6667
• (7 - 8.6667) = -1.6667
• (10 - 8.6667) = 1.3333
• (12 - 8.6667) = 3.3333
• (15 - 8.6667) = 6.3333

Squared deviations:

32.1115 + 13.4445 + 2.7779 + 1.7777 + 11.1109 + 40.1107 = 101.3332

s^2_X = \frac{101.3332}{5} = 20.2666

For Y: Deviations: 3.2833, 2.0833, 3.4833, -1.6167, -3.2167, -4.0167

s^2_Y = \frac{56.3483}{5} = 11.2697

# Verification
var(magnitude)

[1] 20.26667

var(gallagher)

[1] 11.26967

Covariance and Correlation

Covariance: s_{XY} = \frac{\sum(X_i - \bar{X})(Y_i - \bar{Y})}{n-1}

Products of deviations:

• (-5.6667 × 3.2833) = -18.6057
• (-3.6667 × 2.0833) = -7.6387
• (-1.6667 × 3.4833) = -5.8056
• (1.3333 × -1.6167) = -2.1556
• (3.3333 × -3.2167) = -10.7223
• (6.3333 × -4.0167) = -25.4391

Sum = -70.3670

s_{XY} = \frac{-70.3670}{5} = -14.0734

# Verification
cov(magnitude, gallagher)

[1] -14.07333

Correlation: r_{XY} = \frac{s_{XY}}{\sqrt{s^2_X}\sqrt{s^2_Y}} = \frac{-14.0734}{\sqrt{20.2666}\sqrt{11.2697}} = -0.9279

# Verification
cor(magnitude, gallagher)

[1] -0.9312157

OLS Regression (\hat{Y} = \hat{\beta_0} + \hat{\beta_1}X)

Slope coefficient: \hat{\beta_1} = \frac{s_{XY}}{s^2_X} = \frac{-14.0734}{20.2666} = -0.6944

Intercept: \hat{\beta_0} = \bar{Y} - \hat{\beta_1}\bar{X} Steps:

1. -0.6944 × 8.6667 = -6.0181
2. \hat{\beta_0} = 5.7167 - (-6.0181) = 11.7348

# Verification
coef(model)

(Intercept)   magnitude 
 11.7348684  -0.6944079 

R-squared Calculation

Step 1: Calculate predicted values (\hat{Y}):

\hat{Y} = 11.7348 - 0.6944X

# Predicted values
y_hat <- 11.7348 - 0.6944 * magnitude
data.frame(
  Magnitude = magnitude,
  Gallagher = gallagher,
  Predicted = y_hat,
  row.names = 1:6
)

  Magnitude Gallagher Predicted
1         3       9.0    9.6516
2         5       7.8    8.2628
3         7       9.2    6.8740
4        10       4.1    4.7908
5        12       2.5    3.4020
6        15       1.7    1.3188

Step 2: Sum of Squares SST = \sum(Y_i - \bar{Y})^2 = 56.3483 SSR = \sum(\hat{Y}_i - \bar{Y})^2 = 48.5271 SSE = \sum(Y_i - \hat{Y}_i)^2 = 7.8212

R-squared: R^2 = \frac{SSR}{SST} = \frac{48.5271}{56.3483} = 0.8612

# Verification
summary(model)$r.squared

[1] 0.8671626

Visualization - True vs. Estimated Parameters

• True DGP: Y = 12 - 0.8X + ε
• Estimated Model: Y = 11.7348 - 0.6944X

library(ggplot2)

# Create data frame with original data
df <- data.frame(
  magnitude = magnitude,
  gallagher = gallagher
)

# Create sequence for smooth lines
x_seq <- seq(min(magnitude), max(magnitude), length.out = 100)

# Calculate predicted values for both lines
true_dgp <- 12 - 0.8 * x_seq
estimated <- 11.7348 - 0.6944 * x_seq

# Combine into a data frame for plotting
lines_df <- data.frame(
  magnitude = rep(x_seq, 2),
  value = c(true_dgp, estimated),
  Model = rep(c("True DGP", "Estimated"), each = length(x_seq))
)

# Create plot
ggplot() +
  geom_line(data = lines_df, 
            aes(x = magnitude, y = value, color = Model, linetype = Model),
            size = 1) +
  geom_point(data = df, 
             aes(x = magnitude, y = gallagher),
             color = "black", 
             size = 3) +
  scale_color_manual(values = c("red", "blue")) +
  scale_linetype_manual(values = c("dashed", "solid")) +
  labs(
    title = "True DGP vs. Estimated Regression Line",
    subtitle = "Black points show observed data with random noise",
    x = "District Magnitude",
    y = "Gallagher Index",
    caption = "True DGP: Y = 12 - 0.8X + ε\nEstimated: Y = 11.73 - 0.69X"
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    plot.caption = element_text(hjust = 0)
  )

Observations about Model Fit

1. Slope Comparison
   • True slope: -0.8
   • Estimated slope: -0.69
   • The estimated slope is reasonably close to the true parameter
2. Intercept Comparison
   • True intercept: 12
   • Estimated intercept: 11.73
   • The estimated intercept very closely approximates the true value
3. Visual Patterns
   • The lines are nearly parallel, showing good slope recovery
   • Points scatter around both lines due to the random error term (ε)
   • The small sample size (n=6) leads to some imprecision in estimation
   • The estimated line (blue) provides a good approximation of the true DGP (red dashed)
4. Impact of Random Error
   • The scatter of points around the true DGP line reflects the N(0,1) error term
   • This noise leads to the slight differences in estimated parameters
   • With a larger sample, we would expect even closer convergence to true parameters

Interpretation

1. Strong negative correlation (r = -0.93) between district magnitude and electoral disproportionality
2. For each unit increase in district magnitude, the Gallagher index decreases by approximately 0.69 points
3. The model explains 86.12% of the variance in disproportionality
4. The relationship appears strongly linear with moderate scatter
5. The intercept (11.73) represents the expected disproportionality in a hypothetical single-member district system

Study Context

• Data represents simulated observations from a DGP with moderate noise
• Sample shows how increasing district magnitude tends to reduce disproportionality
• Random component reflects other institutional and political factors affecting disproportionality

Limitations

1. Small sample size (n=6)
2. Simulated rather than real-world data
3. Assumes linear relationship
4. Does not account for other institutional features