options(repr.plot.width = 10, repr.plot.height = 10)

Statystyka i Analiza danych¶

Laboratorium 8 - Korelacja i regresja (część 1)¶

Ćwiczenie 1: Wyznaczanie korelacji¶

W ćwiczeniu tym sprawdzimy korelację między dwoma sygnałami. W szczególności sprawdzimy, jak zmiany w definicji oraz we wzajemnym położeniu sygnałów wpływają na wartość korelacji.

Zaczniemy od przygotowania danych. Wykonaj poniższy kod, aby przygotować dane do ćwiczenia.

t <- seq(0,4,0.2)
n <- length(t)

prepareData <- function(x1, x2){
  return(data.frame(t, x1, x2))
}
# Tworzymy różne pary sygnałów
sygnaly_1 <- prepareData(x <- runif(n),1.5*x+runif(n)/2)
sygnaly_2 <- prepareData(sin(t), cos(t))
sygnaly_3 <- prepareData(runif(n), runif(n))
sygnaly_4 <- prepareData(runif(n), 10*runif(n)-5)
sygnaly_5 <- prepareData(x <- runif(n), c(x[n],x[1:n-1]))
sygnaly_6 <- prepareData(sin(t), 5*sin(t))
sygnaly_7 <- prepareData(sin(t)+rnorm(n), sin(t)+rnorm(n))

# Usuwamy niepotrzebne zmienne (w tym funkcje)
rm(n, t, x, prepareData)

Każda ze stworzonych wcześniej ramek danych zawiera informacje o wartościach dwóch sygnałów, x1 i x2 w czasie t. Wyświetl zawartość ramki sygnaly_1.

sygnaly_1

Przypisz wektory sygnaly_1$t, sygnaly_1$x1 i sygnaly_1$x2 do zmiennych t, x1 i x2.

t <- sygnaly_1$t
x1 <- sygnaly_1$x1
x2 <- sygnaly_1$x2

Korzystając z funkcji plot wykonaj wykres rozrzutu (scattered plot), nazwij osie X1 i X2 i nadaj wykresowi tytuł Sygnały 1.

plot(x2 ~ x1, main="Sygnały 1", xlab = "X1", ylab = "X2")

Powstaje pytanie czy wartości sygnałów są ze sobą powiązane - innymi słowy czy są skorelowane? Aby to sprawdzić możemy policzyć pewien punkt odniesienia (wartość średnią) i dla każdego t sprawdzić, po której stronie wartości średniej znajduje się wartość sygnału. Jeżeli dla dwóch sygnałów znak będzie taki sam, to ich iloczyn będzie dodatni, w przeciwnym przypadku ujemny.

Wykonaj poniższy kod, aby zdefiniować funkcję plotWithMeans pozwalającą na umieszczanie na wykresie informacji o przebiegach (w czasie) i wartościach średnich sygnałów.

plotWithMeans <- function(t, x1, x2){
  min <- min(x1, x2)
  max <- max(x1,x2)
  plot(x1 ~ t, ylim=c(min, max), type="o", col="red", ylab="value")
  abline(h=mean(x1), col="red")
  lines(x2 ~ t, col="green",type="o")
  abline(h=mean(x2), col="green")
  legend('topright', c("x1", "x2", "E(x1)", "E(x2)"), lty=1, col=c('red', 'green', 'red',' green'), bty='n', cex=.75)
}

Korzystając z funkcji plotWithMeans stwórz wykres wartości sygnału w zależności od czasu z naniesionymi wartościami średnimi sygnałów.

plotWithMeans(t, x1, x2)

Oblicz wartości średnie sygnałów x1 i x2.

ex1 <- mean(x1)
ex2 <- mean(x2)

Oblicz iloczyny odchyleń wartości sygnałów od średnich (x1-ex1)*(x2-ex2).

prod <- (x1-ex1)*(x2-ex2)

Uśredniając takie iloczyny otrzymujemy kowariancję dwóch sygnałów. Oblicz kowariancję (z próby) dla x1 i x2. Zweryfikuj poprawność obliczeń korzystając z wbudowanej funkcji cov.

kow <- sum(prod)/(length(prod)-1)
kow
cov(x1, x2)

Kowariancja jest zależna od wartości bezwględnych. Jej podzielenie przez iloczyn odchyleń standardowych normuje wartość współczynnika do przedziału od -1 do 1 - jest to korelacja. Oblicz wartość współczynnika korelacji $\frac{cov(x1, x2)}{sd(x1)\cdot sd(x2)}$.

kor <- kow/(sd(x1)*sd(x2))
kor
cor(x1, x2)

Na podstawie wykonanych wcześniej poleceń uzupełnij funkcję analyzeCorr, która przyjmie jako argumenty t, x1, x2 i tytuł (title), stworzy wykres rozrzutu i wykres zależności sygnałów od czasu wraz z średnimi (wykresy powinny być wyświetlone obok siebie) i zwróci wartość współczynnika korelacji dla x1 i x2. Do wykreślenia wykresów obok siebie możesz skorzystać z polecenia par(mfrow=c(1,2)).

analyzeCor <- function(t, x1, x2, title){
    # Ustawiamy dwa wykresy obok siebie
    prev  <- par(mfrow=c(1,2))
    # Wyznaczamy kowariancję i korelację
    ex1 <- mean(x1)
    ex2 <- mean(x2)
    prod <- (x1-ex1)*(x2-ex2)
    cov <- sum(prod)/(length(prod)-1)
    cor <- cov/(sd(x1)*sd(x2))
    # 1. wykres rozrzutu 
    plot(x2 ~ x1, xlab="X1", ylab="X2", main=paste("r =",round(cor, 2)))
    # 2. wkres plotWithMeans
    plotWithMeans(t, x1, x2)
    # Przywracamy poprzednie ustawienia
    par(prev)
    return(cor)
}

Sprawdzić działanie funkcji analyzeCorr.

kor  <- analyzeCor(t, x1, x2, "sygnały 1")
kor

Za pomocą funkcji analyzeCorr przeanalizuj wszystkie pary sygnałów wykorzystując poniższą pętlę.

cor  <- vector(length = 7)
names(cor)  <- paste("sygnaly", 1:7, sep="_")
for (name in names(cor)){
  signals <- eval(as.name(name))
  cor[name] <- analyzeCor(signals$t, signals$x1, signals$x2, name)
}

A następnie wyświetl współczynniki korelacji dla poszczególnych par sygnałów (zmnienna cor).

cor

Porównaj powyższe wyniki z wartościami uzykanymi za pomocą funkcji wbudowanej cor.

cor  <- vector(length = 7)
names(cor)  <- paste("sygnaly", 1:7, sep="_")
for (name in names(cor)){
  signals <- eval(as.name(name))
  cor[name] <- cor(signals$x1, signals$x2)
}
cor

Ćwiczenie 2: Wyznaczanie współczynników modelu regresji¶

W ćwiczeniu tym wyznaczymy współczynniki w modelu regresji stosując metodę najmnieszych kwadratów.

Podobnie jak poprzednio, zaczniemy od przygotowania danych. Wykonaj poniższy kod, aby przygotować dane do ćwiczenia.

regresja <- read.csv(url("http://www.cs.put.poznan.pl/swilk/siad/8-cw2.csv"), sep=";")

Ramka zawiera trzy kolumny, x, yp=f(x) i y=f(x)+epsilon. Postaramy się wyestymować parametry a i b, tak aby jak najlepiej przybliżyć wartość y. Pierwsze podejście będzie polegać na ręcznym dopasowaniu wartości a^ i b^, tak, aby suma różnic kwadratów była jak najmniejsza.

Uzupełnij definicję funkcji plotReg, która wykreśla yp, y oraz y^=a^*x+b^ i zwraca wartość sumy kwadratów różnic (y^-y)^2.

plotReg <- function(reg, at, bt){
    # Oblicz y^=a^*x+b^:
    yt = reg$x*at + bt
    # Oblicz sumę (Y^-Y)^2:
    ss  <-  sum((yt-reg$y)^2)
  
    plot(reg$yp ~ reg$x, t="l", col="red", xlab="x", ylab="")
    points(reg$y ~ reg$x, col="blue")
    lines(yt ~ reg$x, t="l", col="green")
  
    # Funkcja powinna zwracać sumę kwadratóW różnic:
    return(ss)
}

Spróbuj dopasować wartości a^ i b^ korzystając ze zdefiniowanej powyżej funkcji.

plotReg(regresja, 3, 2.9)

Następnie wykorzystamy funkcję lm. Pierwszym argumentem tej metody jest formula, którą zdefinijemy jako y~x, drugim jest data=regresja. Otrzymany model regresji zapisz do zmiennej model.

model <- lm(y ~ x, regresja)

Korzystając ze zmiennej model odczytaj 'a^' i 'b^'.

model

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

Coefficients:
(Intercept)            x  
      2.913        3.001

Obszerniejsze wyniki można uzyskać korzystając z funkcji summary.

summary(model)

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

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4594 -0.8644 -0.2456  0.7622  2.2902 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.91253    0.22896   12.72 1.14e-12 ***
x            3.00143    0.03873   77.50  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.032 on 26 degrees of freedom
Multiple R-squared:  0.9957,	Adjusted R-squared:  0.9955 
F-statistic:  6006 on 1 and 26 DF,  p-value: < 2.2e-16

Ćwiczenie 3: Testowanie istotności statystycznej korelacji¶

W tym ćwiczeniu sprawdzimy, czy korelacja między zmiennymi jest statystycznie istotna. W tym celu skorzystamy ze statystyki t.

Zwyczajowo już zaczynamy od załadowania danych wykonyjąc poniższy kod.

dzieci <- read.csv(url("http://www.cs.put.poznan.pl/swilk/siad/8-cw3.csv"), sep=";")

Dane te zawierają informację o wieku w latach (x) i wzroście w cm (y) dla 15-osobowej grupy wylosowanej z populacji dzieci i młodzieży. Wyświetl zawartości ramki dzieci.

dzieci

Korzystajac z funkcji plot utwórz wykres rozrzutu dla tych danych. Opisz odpowiednio osie wykresu.

plot(dzieci$y ~ dzieci$x, xlab="wiek [lata]", ylab="wzrost [cm]")

Oblicz średnie i odchylenia standardowe obu zmiennych oraz liczbę obserwacji.

ex <- mean(dzieci$x)
ey <- mean(dzieci$y)
sdx <- sd(dzieci$x)
sdy <- sd(dzieci$y)
n <- nrow(dzieci)

Oblicz samodzielnie współczynnik korelacji korzystając ze wzoru $\frac{cov(x, y)}{sd(x) \cdot sd(y)}$.

sum((dzieci$x - ex)*(dzieci$y - ey))/((n - 1)*sdx*sdy)

Oblicz także współczynnik korelacji korzystając z funkcji cor.

cor  <- cor(dzieci$x, dzieci$y)
cor

Aby zweryfikować hipotezę mówiącą, że współczynnik korelacji jest niezerowy (test dwustronny), oblicz wartość statystyki $t=\frac{r}{\sqrt{1-r^2}}\cdot \sqrt{n-2}$, gdzie $r$ to współczynnik korelacji dla próby.

t <- cor/sqrt(1-cor^2)*sqrt(n-2)
t

Oblicz p-wartość dla wyznaczonej statystyki pamiętając, że mamy do czynienia z testem dwustronnym. Jaki będzie wynik testowania hipotez (zakładając poziom istotności równy 0.05)?

2*pt(t, lower.tail=F, df=n-2)
cor.test(dzieci$x, dzieci$y)

	Pearson's product-moment correlation

data:  dzieci$x and dzieci$y
t = 45.263, df = 13, p-value = 1.086e-15
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9902420 0.9989805
sample estimates:
      cor 
0.9968424

Wyznacz model regresji korzystając z funkcji lm.

model <- lm(y ~ x, data=dzieci)
model

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

Coefficients:
(Intercept)            x  
     88.689        4.305

Utwórz wykres modelu regresji korzystając z funkcji abline (wcześniej narysuj ponownie wykres rozrzutu). W jakich warunkach uzyskany model może mieć zastosowanie?

plot(y ~ x, data=dzieci, xlab="wiek [lata]", ylab="wzrost [cm]")
abline(model, col="red")

t	x1	x2
<dbl>	<dbl>	<dbl>
0.0	0.3902706	0.7796120
0.2	0.3424317	0.6154584
0.4	0.4504058	0.7495732
0.6	0.9970429	1.7636664
0.8	0.6309133	1.2215927
1.0	0.8776847	1.5925431
1.2	0.4972250	1.0926869
1.4	0.8386165	1.5530709
1.6	0.9407478	1.4282569
1.8	0.8383732	1.7430421
2.0	0.9679630	1.5709541
2.2	0.9501762	1.4372536
2.4	0.1781327	0.2727217
2.6	0.8329675	1.3631238
2.8	0.5286438	1.1172406
3.0	0.5700783	0.8823672
3.2	0.1364164	0.2597273
3.4	0.2287090	0.5438605
3.6	0.8686166	1.6299434
3.8	0.9526149	1.8567394
4.0	0.4706623	0.9279454

x	y
<dbl>	<int>
7.0	120
8.0	122
9.0	125
10.0	131
11.0	135
11.5	140
12.0	142
13.0	145
14.0	150
15.0	154
16.0	159
17.0	162
18.0	164
18.5	168
19.0	170