Statystyka i Analiza danych

Laboratorium 9 - Korelacja i regresja (część 2)

Ćwiczenie 1: Regresja liniowa w R

Celem ćwiczenia jest zapoznanie się z funkcją lm() wykonującą regresję liniową w R.

Zacznijmy od wczytania poniższej funkcji generującej prostą ramkę danych:

In [1]:
generujDane <- function(n, a, b, sigma) {
  x <- runif(n)
  y <- a * x + b + sigma * rnorm(n)
  plot(y ~ x, t="p", xlab="X", ylab="Y", main="Regression Data")
  abline(b, a, lty="dashed")
  data.frame(X=x, Y=y)
}

Użyj powyższej funkcji aby wygenerować zbiór danych Regresja składający się z 50 obserwacji z prawdziwymi współczynnikami a=3, b=-2 i szumem sigma=1.

In [2]:
Regresja <- generujDane(50, 3, -2, 1)

Wyznacz współczynniki regresji liniowej, wyświetl zmienną model i porównaj współczynniki z prawdziwymi wartościami a, b.

In [3]:
model <- lm(Y ~ X, Regresja)
model

Call:
lm(formula = Y ~ X, data = Regresja)

Coefficients:
(Intercept)            X  
     -1.874        2.701

Przypisz współczynniki do zmiennej coeffs i wyznacz wektor $\hat{Y}$ (ozn. Yp) wartości zmiennej objaśnianej wyznaczone z modelu liniowego.

In [4]:
coeffs <- model$coefficients
Yp <- Regresja$X * coeffs[2] + coeffs[1]

Wyznacz i wyświetl podsumowanie modelu. Zinterpretuj wypisywane wartości

In [5]:
modelSummary = summary(model)
modelSummary

Call:
lm(formula = Y ~ X, data = Regresja)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.1547 -0.6561 -0.1227  0.6611  2.4383 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.8744     0.2829  -6.625 2.76e-08 ***
X             2.7010     0.5192   5.203 4.03e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.016 on 48 degrees of freedom
Multiple R-squared:  0.3606,	Adjusted R-squared:  0.3473 
F-statistic: 27.07 on 1 and 48 DF,  p-value: 4.03e-06

Do poszczególnych elementów podsumowania można uzyskać dostęp poprzez funkcje:

  • współczynniki modelu z błędami standardowymi, wartościami statystyki t i jej p-wartościami:
In [6]:
modelSummary$coefficients
A matrix: 2 × 4 of type dbl
EstimateStd. Errort valuePr(>|t|)
(Intercept)-1.8743620.2829083-6.6253352.762025e-08
X 2.7009690.5191560 5.2026164.029684e-06
  • Współczynnik determinacji $R^2$:
In [7]:
modelSummary$r.squared
0.360573037091662
  • Statystyka $F$:
In [8]:
F <- modelSummary$fstatistic
F
value
27.0672129646819
numdf
1
dendf
48

Sprawdź, czy model jest istotny statystycznie za pomocą statystyki F i porównaj z p-wartością zwracaną przez modelSummary:

In [9]:
1 - pf(F[1], F[2], F[3])
value: 4.02968391854941e-06

Wyznacz również SST, SSE i SSR:

In [10]:
yMean <- mean(Regresja$Y)
SST <- sum((Regresja$Y - yMean)^2)
SSR <- sum((Yp - yMean)^2)
SSE <- sum((Regresja$Y - Yp)^2)
SST - SSR - SSE
1.4210854715202e-14

Na koniec możesz wyświetlić wykresy diagnostyczne:

In [11]:
plot(model)

Czyszczenie przestrzeni roboczej (usunięcie wszystkich zmiennych, funkcji, itp.)

In [12]:
rm(list = ls())

Ćwiczenie 2: Dane medyczne

Wysunięto hipotezę, że istnieje związek pomiędzy czasem działania pewnego leku u chorych na zaburzenia układu krążenia a aktywnością pewnego enzymu. Losowa próba dała następujące wyniki (patrz tabela). Dla zmiennych x i y należy:

  • wyliczyć współczynnik korelacji Pearsona
  • utworzyć wykres rozrzutu
  • obliczyć funkcję regresji i dodaj do wykresu linię trendu
  • Czy odrzucenie danych o pewnych obserwacjach mogłoby poprawić wyniki?
  • Czy jeszcze jakaś transformacja zbioru danych pozwoliłaby na dalsze polepszenie wyników ?

Zaczynamy od załadowania danych wykonując poniższy kod:

In [13]:
Dane <- read.csv(url("http://www.cs.put.poznan.pl/wkotlowski/siad/9-cw2.csv"), sep=";")

Dane te zawierają informację o płci, czasie działania leku (x) i aktywności enzymu (y). Wyświetl tabelę z danymi:

In [14]:
Dane
A data.frame: 22 × 3
plecxy
<fct><dbl><dbl>
k15.644.5
k10.644.0
k16.925.8
k15.039.3
k 6.721.8
k12.037.3
k16.043.5
k13.021.3
k 9.732.8
k 7.843.0
k15.0 7.1
k14.0 5.0
m21.723.1
m14.529.6
m17.327.1
m22.423.6
m20.925.2
m20.324.2
m24.722.2
m15.328.8
m18.626.9
m23.622.5

Wyznacz współczynnik korelacji Pearsona:

In [15]:
wspKor <- cor(Dane$x, Dane$y)
wspKor
-0.302252739234036

Korzystajac z funkcji plot utwórz wykres rozrzutu dla tych danych. Opisz odpowiednio osie wykresu. Wyznacz funkcję regresji i dodaj ją do wykresu:

In [16]:
plot(Dane$y ~ Dane$x, main = "Dane medyczne", xlab = "x (aktywność enzymu)", ylab = "y (czas działania leku)")
model <- lm(y ~ x, Dane)
abline(model)

Sprawdź, czy wynik jest istotny na poziomie istotności alfa=0.05 za pomocą funkcji summary lub ręcznie wyznaczając p-wartość statystyki F:

In [17]:
model.summary = summary(model)
F <- model.summary$fstatistic
pValue <- 1 - pf(F[1], F[2], F[3])
model.summary
pValue

Call:
lm(formula = y ~ x, data = Dane)

Residuals:
     Min       1Q   Median       3Q      Max 
-24.4163  -1.2270   0.0402   5.0721  16.1317 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  38.5863     7.7131   5.003 6.83e-05 ***
x            -0.6550     0.4619  -1.418    0.172    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.48 on 20 degrees of freedom
Multiple R-squared:  0.09136,	Adjusted R-squared:  0.04592 
F-statistic: 2.011 on 1 and 20 DF,  p-value: 0.1716
value: 0.171572744235917

Czy odrzucenie części danych mogłoby poprawić wyniki?

In [18]:
DaneK <- Dane[Dane$plec == "k",]
plot(DaneK$y ~ DaneK$x, main = "Dane medyczne (kobiety)", xlab = "x (aktywność enzymu)", ylab = "y (czas działania leku)")
model <- lm(y ~ x, DaneK)
abline(model)
summary(model)

Call:
lm(formula = y ~ x, data = DaneK)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.868  -9.588   3.781  10.936  15.343 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  36.0944    17.4449   2.069   0.0654 .
x            -0.4447     1.3327  -0.334   0.7455  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 14.79 on 10 degrees of freedom
Multiple R-squared:  0.01101,	Adjusted R-squared:  -0.08789 
F-statistic: 0.1114 on 1 and 10 DF,  p-value: 0.7455
In [19]:
DaneM <- Dane[Dane$plec == "m",]
plot(DaneM$y ~ DaneM$x, main = "Dane medyczne (mezczyzni)", xlab = "x (aktywność enzymu)", ylab = "y (czas działania leku)")
model <- lm(y ~ x, DaneM)
abline(model)
summary(model)

Call:
lm(formula = y ~ x, data = DaneM)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.88466 -0.16592  0.06522  0.40484  0.61180 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 40.35580    1.10106   36.65 3.37e-10 ***
x           -0.75443    0.05452  -13.84 7.19e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5617 on 8 degrees of freedom
Multiple R-squared:  0.9599,	Adjusted R-squared:  0.9549 
F-statistic: 191.5 on 1 and 8 DF,  p-value: 7.192e-07

Czyścimy przestrzeń roboczą:

In [20]:
rm(list = ls())

Ćwiczenie 3: Szacowanie czasu pracy programistów

Baza danych historycznych obejmuje oszacowane przez programistów czasy pracy oraz faktyczny rozmiar programu. Sprawdź, czy te dwie zmienne są zależne i przeprowadź analizę regresji. Powyższe zadanie ilustruje jeden z elementów metody PSP.

Wczytaj i wyświetl dane z pliku:

In [21]:
Programy <- read.csv(url("http://www.cs.put.poznan.pl/wkotlowski/siad/9-cw3.csv"), sep=";")
Programy
A data.frame: 10 × 2
rozmiarczas
<int><int>
186130
699650
132 99
272150
291128
331302
199 95
1890945
788368
1601961

Teraz dołączymy Programy do ścieżki przeszukiwania nazw zmiennych. Dzięki temu możemy odwoływać się do zmiennych w Programy bez konieczności każdorazowego zaznaczania ramki danych

In [22]:
attach(Programy)

Np. wyznaczenie współczynnika korelacji Pearsona wygląda teraz tak:

In [23]:
wspKor <- cor(rozmiar, czas)
wspKor
0.954496574104683

Utwórz wykres rozrzutu (czas w funkcji rozmiaru), oblicz funkcję regresji i dodaj ją do wykresu:

In [24]:
plot(czas ~ rozmiar, main = "Programy", xlab = "rozmiar", ylab = "czas")
model <- lm(czas ~ rozmiar)
abline(model)

Sprawdź, czy wynik jest istotny na poziomie istotności alfa=0.05 za pomocą funkcji summary

In [25]:
model.summary = summary(model)
F <- model.summary$fstatistic
pValue <- 1 - pf(F[1], F[2], F[3])
pValue
value: 1.77517178131525e-05

Zinterpretuj współczynniki modelu

In [26]:
model

Call:
lm(formula = czas ~ rozmiar)

Coefficients:
(Intercept)      rozmiar  
    45.9358       0.5273

Odłącz ramkę danych od ścieżki i wyczyść przestrzeń roboczą:

In [27]:
detach(Programy)
rm(list = ls())

Ćwiczenie 4: Analiza regresji wielu zmiennych (wieloraka)

Zadanie demonstracyjne, które pokaże przykład analizy regresji wielu zmiennych i różnicę między testem t na istotoność zmiennych w modelu, a globalnym testem F.

Wczytaj dane z pliku 9-cw4.csv i wyświetl tabelę danych. Jest to zbiór danych ekonomicznych, gdzie Exports jest zmienną objaśnianą (Y), a pozostałe są zmiennymi objaśniającymi (X)

In [28]:
Ekonom <- read.csv(url("http://www.cs.put.poznan.pl/wkotlowski/siad/9-cw4.csv"), sep=";")
Ekonom
A data.frame: 67 × 5
ExportsM1LendPriceExchange
<dbl><dbl><dbl><dbl><dbl>
2.65.1 7.81142.2
2.64.9 8.01162.2
2.75.1 8.11172.2
3.05.1 8.11222.2
2.95.1 8.11242.2
3.15.2 8.11282.2
3.25.1 8.31322.1
3.75.2 8.81332.2
3.65.3 8.91332.2
3.45.4 9.11342.2
3.75.7 9.21352.2
3.65.7 9.51362.2
4.15.910.31402.2
3.55.810.61472.2
4.25.711.31502.2
4.35.812.11512.2
4.26.012.01512.2
4.16.011.41512.1
4.66.011.11532.1
4.46.011.01542.1
4.56.111.31542.1
4.66.012.61542.1
4.66.113.61552.1
4.26.713.61552.1
5.56.214.31562.1
3.76.314.31562.1
4.97.013.71592.1
5.27.012.71612.1
4.96.612.61612.2
4.66.413.41612.1
⋮⋮⋮⋮⋮
4.87.110.11542.1
5.47.010.01542.1
5.07.510.21542.1
5.27.411.01532.0
4.77.411.01522.1
5.17.310.71522.2
4.97.610.21522.2
4.97.810.01512.2
5.37.8 9.81522.2
4.88.2 9.31522.2
4.98.2 9.31522.2
5.18.3 9.51522.1
4.38.3 9.21502.1
4.98.0 9.11472.1
5.38.2 9.01472.1
4.88.2 9.01462.1
5.38.0 8.91452.1
5.08.1 9.01452.1
5.18.1 9.01462.1
4.88.1 9.01472.1
4.88.1 8.91472.1
5.28.6 8.91472.1
4.98.8 9.01462.1
5.58.4 9.11472.1
4.38.2 9.01462.1
5.28.3 9.21462.1
4.78.3 9.61462.1
5.48.410.01462.1
5.28.310.01472.1
5.68.210.11462.2

Utwórz model liniowy zmiennej Exports w funkcji pozostałych zmiennych:

In [29]:
model <- lm(Exports ~ ., Ekonom) # skrócony zapis dla lm(Exports ~ M1 + Lend + Price + Exchange, Ekonom)
model

Call:
lm(formula = Exports ~ ., data = Ekonom)

Coefficients:
(Intercept)           M1         Lend        Price     Exchange  
  -4.047103     0.368372     0.002669     0.036862     0.268003

Wyświetl podsumowanie modelu:

In [30]:
modelSummary <- summary(model)
modelSummary

Call:
lm(formula = Exports ~ ., data = Ekonom)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.92502 -0.12779 -0.03198  0.18502  0.91181 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4.047103   2.329192  -1.738   0.0873 .  
M1           0.368372   0.062276   5.915 1.55e-07 ***
Lend         0.002669   0.047919   0.056   0.9558    
Price        0.036862   0.009324   3.954   0.0002 ***
Exchange     0.268003   0.928517   0.289   0.7738    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3357 on 62 degrees of freedom
Multiple R-squared:  0.8251,	Adjusted R-squared:  0.8138 
F-statistic:  73.1 on 4 and 62 DF,  p-value: < 2.2e-16

O ile sam model jest istotny statystycznie na poziomie alfa=0.05 (p-wartość statystyki F < 2.2e-16), o tyle niektóre ze zmiennych nie są istotne statystycznie (p-wartość statystyki t < 0.05). Poniższa analiza może nam posłużyć do selekcji zmiennych istotnych statystycznie.

Odrzucamy więc zmienne nieistotne statystycznie i pozostawiamy tylko Exports, M1 i Price

In [31]:
Ekonom2 <- Ekonom[, c("Exports", "M1", "Price")]

Ponownie wykonujemy analizę regresji:

In [32]:
model <- lm(Exports ~ ., Ekonom2) 
model
summary(model)

Call:
lm(formula = Exports ~ ., data = Ekonom2)

Coefficients:
(Intercept)           M1        Price  
   -3.42296      0.36142      0.03703  

Call:
lm(formula = Exports ~ ., data = Ekonom2)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.93106 -0.12490 -0.02465  0.18140  0.90508 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.422957   0.540853  -6.329 2.75e-08 ***
M1           0.361417   0.039246   9.209 2.45e-13 ***
Price        0.037033   0.004094   9.046 4.70e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3306 on 64 degrees of freedom
Multiple R-squared:  0.8248,	Adjusted R-squared:  0.8193 
F-statistic: 150.7 on 2 and 64 DF,  p-value: < 2.2e-16