Kowariancja

Analiza kowariancji to często stosowana metoda statystyczna łącząca w sobie elementy analizy wariancji, korelacji i regresji.

Główny cel metody jest podobny do analizy wariancji: dać odpowiedź czy analizowany czynnik doświadczalny wpływa w sposób istotny na badaną cechę. Zatem metodę tę stosujemy gdy na wybraną cechę działa czynnik doświadczalny na kilku poziomach. Dodatkowo, zaletą metody analizy kowariancji jest możliwość wyeliminowania wpływu innej cechy mającej wpływ na cechę badaną.
Tą inną cechę nazywać będziemy zmienną towarzyszącą.

1. Nasze dane (tym razem nasze dane nie są czytane jak do tej pory z pliku o rozszerzeniu .csv, ale tworzymy je jako wektory , to c przed nawiasami oznacza, ze dane są danymi wektorowymi)

paramG<-c(20, 23, 30, 25, 34, 40, 27, 38, 24, 31, 19, 26, 33, 35, 30, 31, 34, 28, 23, 22)

paramD<-c(5, 10, 12, 9, 23, 21, 14, 18, 6, 13, 7, 12, 27, 24, 18, 22, 26, 21, 14, 9)

2. Nie wiem co to pozwalająca na powtórzenie wektora określoną ilość razy, np. rep(c(1,2),2) jest równoważne c(1,2,1,2). Faktor (factor) jest specjalną strukturą, przechowującą (oprócz szeregu danych) informacje o powtórzeniach takich samych wartości oraz o zbiorze unikalnych wartości tego ciągu.

Dodajemy nową daną „lek” (która również jest wektorem) i przypisujemy jej następujące wartości A i B, przy czym obie te wartości mają być powtórzone po 10 razy.

lek<-as.factor(c(rep("A",10),rep("B",10)))

3. Wyświetlenie danych

paramG paramD lek

1 20 5 A

2 23 10 A

3 30 12 A

4 25 9 A

5 34 23 A

6 40 21 A

7 27 14 A

8 38 18 A

9 24 6 A

10 31 13 A

11 19 7 B

12 26 12 B

13 33 27 B

14 35 24 B

15 30 18 B

16 31 22 B

17 34 26 B

18 28 21 B

19 23 14 B

20 22 9 B

4. Graficzne przedstawienie poszczególnych wartości parametru G zarówno dla leku A , jak i dla leku B

stripchart(paramG~lek,vertical=TRUE,method="stack")

Wykres przedstawia pionowo wartości parametru G w zależności od stosowanego leku (lek A i lek B).

Widoczne jest zróżnicowanie wartości parametru G w przypadku stosowanych leków.

Wartości średnie parametru G w poszczególnych grupach

5. Funkcja tapply (Grupuje ona dane liczbowe w grupy odpowiadające poszczególnym danym kategorycznym, a następnie przeprowadza na tych grupach określone operacje).

tapply(paramG,lek,mean)

A B

29.2 28.1

6. Średnie w każdej grupach testowanych leków

by(paramG,lek,mean)

lek: A

[1] 29.2

------------------------------------------------------------

lek: B

[1] 28.1

MOŻNA TU POLICZYĆ RESZTY I SPRAWDZIĆ ICH ROZKŁAD NORMALNOŚCI

7. Shapiro test dla parametru G w grupach testowanych leków ( dla leku A i dla leku B) Wysoka wartość p świadczy o tym, że nie ma podstaw do odrzucenia hipotezy o normalności rozkładu tych danych.

H₀: Występuje normalność rozkładu danych.

H₁: Dane nie podlegają rozkładowi normalnemu.

by(paramG,lek,shapiro.test)

lek: A

Shapiro-Wilk normality test

data: dd[x, ]

W = 0.9571, p-value = 0.7526

------------------------------------------------------------

lek: B

Shapiro-Wilk normality test

data: dd[x, ]

W = 0.9478, p-value = 0.6429

Zarówno w przypadku leku A, jak i leku B, p- wartości są większe od α= 0,05. Nie ma zatem podstaw do odrzucenia hipotezy zerowej, co oznacza, że wśród danych występuje normalność rozkładu.

8. Odchylenie standardowe

by(paramG,lek,sd)

lek: A

[1] 6.613118

------------------------------------------------------------

lek: B

[1] 5.466057

9. Analiza jednorodności wariancji sprawdzający, czy wariancje obu prób różnią się istotnie między sobą

H₀: Istnieje różnorodność wariancji.

H₁: Wariancję różnią się między sobą.

var.test(paramG~lek)

F test to compare two variances

data: paramG by lek

F = 1.4637, num df = 9, denom df = 9, p-value = 0.5794

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

0.3635726 5.8930134

sample estimates:

ratio of variances

1.463741

P- wartość dla przeprowadzonego testu wynosi 0,5794 i jest znacznie większa od α = 0,05. Zatem pomiędzy wariancjami obu prób (wariancje dla parametru G przy stosowaniu leku A
i wariancje dla parametru G przy stosowaniu leku B) nie ma istotnie statystycznie różnic.

Występuje jednorodność wariancji, zatem możemy wykonać test T-Studenta.

10. Ustalenie istotności statystycznej pomiędzy wariancjami średnich parametru G przy stosowaniu leku A i leku

Test t- Studenta pozwala ustalić czy różnica pomiędzy średnimi jest istotna statystycznie, czyli czy możemy odrzucić hipotezę zerową

przy ustalonym poziomie istotności

H₀: Występuje równość wariancji.

H₁: Istnieją różnice pomiędzy wariancjami.

t.test(paramG~lek,var.equal=TRUE)

Two Sample t-test

data: paramG by lek

t = 0.4054, df = 18, p-value = 0.6899

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-4.600089 6.800089

sample estimates:

mean in group A mean in group B

29.2 28.1

Żeby można było przeprowadzić test T-Studenta trzeba wcześniej zrobić Shapiro-Wilka i F-test.

P- wartość przeprowadzonego testu wyniosła 0,6899 i jest ona większa od ustalonego poziomu istotności wynoszącego α =0,05. Oznacza to, że nie ma podstaw do odrzucenia hipotezy zerowej. Zatem pomiędzy wariancjami średnich wartości parametru G przy stosowaniu leku A i leku B nie występują statystycznie istotne różnice.

11. Analiza wariancji

H₀: Wartości parametru G w przypadku stosowaniu obu leków są jednakowe (→leki nie różnią się).

H₁: Wartości parametru G w przypadku stosowaniu obu leków różnią się (→leki różnią się).

summary(aov(paramG~lek))

Df Sum Sq Mean Sq F value Pr(>F)

lek 1 6.0 6.05 0.164 0.69

Residuals 18 662.5 36.81

P- wartość przeprowadzonego testu wyniosła 0,69 i jest ona większa od ustalonego poziomu istotności wynoszącego α =0,05. Oznacza to, że nie ma podstaw do odrzucenia hipotezy zerowej. Zatem nie istnieją różnice pomiędzy lekami, które w sposób statystycznie istotny wpływałyby na wartości parametru G.

12. Wyznaczenie współczynnika korelacji pomiędzy parametrami D i G

cor(paramG,paramD)

[1] 0.802848

Współczynnik korelacji liniowej przyjmuje zawsze wartości w zakresie [-1,1]. Im większa wartość bezwzględna współczynnika, tym większa jest zależność liniowa między zmiennymi.

0→ brak korelacji

1→ silna korelacja (jeśli jeden czynnik rośnie, to drugi też)

-1→ silna korelacja (jeśli czynnik X rośnie to czynnik Y maleje albo odwrotnie)

13. Wyznaczenie istotności współczynnika korelacji pomiędzy parametrami D i G (to u góry nie jest w sumie potrzebne, bo w tej funkcji wyświetla się też sam współczynnik korelacji).

H₀: Parametr G i parametr D nie są od siebie zależne.

H₁: Istnieje współzależność pomiędzy parametrami G i D.

cor.test(paramG,paramD)

Pearson's product-moment correlation

data: paramG and paramD

t = 5.7133, df = 18, p-value = 2.038e-05

alternative hypothesis: true correlation is not equal to 0

95 percent confidence interval:

0.5588867 0.9189035

sample estimates:

cor

0.802848

P- wartość przeprowadzonego testu wyniosła 2,038 * 10^-5 i jest ona mniejsza od ustalonego poziomu istotności wynoszącego α =0,05. Oznacza to, że hipotezę zerowa i przyjmujemy hipotezę alternatywną. Zatem istnieje korelacja pomiędzy wartościami parametru G i D i jest ona statystycznie istotna.

14. Wykres rozrzutu zmiennych empirycznych

plot(paramD,paramG)

Wykres przedstawiający zależności pomiędzy wartościami parametru D i wartościami parametru G.

15. Przypisanie lekom A i B odpowiednich parametrów G i D

paramD_lekA<-paramD[lek=="A"]

paramD_lekB<-paramD[lek=="B"]

paramG_lekA<-paramG[lek=="A"]

paramG_lekB<-paramG[lek=="B"]

16. Graficzne przedstawienie przypisania lekom A i B odpowiednich parametrów

plot(paramD,paramG,type="n")

points(paramD_lekA,paramG_lekA,pch="A")

points(paramD_lekB,paramG_lekB,pch="B",col="red")

Na osi X parametr D, na osi Y parametr G, punkty A odpowiadają lekowi A, punkty B są dla B i dodatkowo B oznaczone kolorem czerwonym. To jest ten sam wykres co wyżej, tylko zamiast punktów są literki.

16. Oszacowanie odrębnych modeli liniowych

• Dla leku A

lm(paramG_lekA~paramD_lekA)

Call:

lm(formula = paramG_lekA ~ paramD_lekA)

Coefficients:

(Intercept) paramD_lekA

16.4226 0.9754

• Dla leku B

lm(paramG_lekB~paramD_lekB)

Call:

lm(formula = paramG_lekB ~ paramD_lekB)

Coefficients:

(Intercept) paramD_lekB

15.1087 0.7217

Analiza kowariancji

17. Sprawdzenie równości współczynników liniowych dla skonstruowanych modeli

H₀: Współczynniki kierunkowe obu modeli liniowych są równe (czyli nie istnieje interakcja pomiędzy parametrem G, a parametrem D i rodzajem zastosowanego leku).

H₁: Występuje różnica pomiędzy kierunkowymi modeli (czyli istnieje interakcja pomiędzy parametrem G, a parametrem D i rodzajem zastosowanego leku).

18. Konstrukcja modelu (zawiera on interakcje między parametrem G,
a parametrem D i rodzajem leku)

model<-lm(paramG~paramD*lek)

Call:

lm(formula = paramG ~ paramD * lek)

Coefficients:

(Intercept) paramD lekB paramD:lekB

16.4226 0.9754 -1.3139 -0.2536

19. Wyznaczenie parametrów modelu i szacowanie ich istotności

• Podsumowanie dla modelu liniowego

summary(model)

Call:

lm(formula = paramG ~ paramD * lek)

Residuals:

Min 1Q Median 3Q Max

-4.8562 -1.7500 0.0696 1.8982 4.0207

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 16.4226 2.0674 7.944 6.08e-07 ***

paramD 0.9754 0.1446 6.747 4.69e-06 ***

lekB -1.3139 3.1310 -0.420 0.680

paramD:lekB -0.2536 0.1893 -1.340 0.199

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.622 on 16 degrees of freedom

Multiple R-squared: 0.8355, Adjusted R-squared: 0.8046

F-statistic: 27.09 on 3 and 16 DF, p-value: 1.655e-06

• Podsumowanie dla modelu liniowego funkcją anova

summary.aov(model)

Df Sum Sq Mean Sq F value Pr(>F)

paramD 1 430.9 430.9 62.689 6.34e-07 ***

lek 1 115.3 115.3 16.774 0.000844 ***

paramD:lek 1 12.3 12.3 1.795 0.199066

Residuals 16 110.0 6.9

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

P- wartość przeprowadzonego testu wyniosła 0,199 i jest ona większa od ustalonego poziomu istotności wynoszącego α =0,05. Oznacza to, że nie ma podstaw do odrzucenia hipotezy zerowej. Zatem nie istnieje interakcja pomiędzy wartością parametru G, a wartością parametru D
i rodzajem zastosowanego leku.

20. Konstrukcja modelu mniej złożonego (bez uwzględnienia interakcji między parametrem G, parametrem D i rodzajem leku)

model2<-lm(paramG~lek+paramD)

model2

Call:

lm(formula = paramG ~ lek + paramD)

Coefficients:

(Intercept) lekB paramD

18.3600 -5.1547 0.8275

21. Wyznaczenie parametrów modelu i szacowanie ich istotności

summary(model2)

Call:

lm(formula = paramG ~ lek + paramD)

Residuals:

Min 1Q Median 3Q Max

-3.6348 -2.5099 -0.2038 1.8871 4.7453

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 18.3600 1.5115 12.147 8.35e-10 ***

lekB -5.1547 1.2876 -4.003 0.000921 ***

paramD 0.8275 0.0955 8.665 1.21e-07 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.682 on 17 degrees of freedom

Multiple R-squared: 0.817, Adjusted R-squared: 0.7955

F-statistic: 37.96 on 2 and 17 DF, p-value: 5.372e-07

P- wartość wyznaczona dla leku i dla parametru D wynoszą odpowiednio 0,000921 i 1,21 * 10^-7, zatem obie są mniejsze od ustalonego poziomu istotności wynoszącego α= 0,05. Oznacza to, że należy odrzucić hipotezę zerową i jako właściwą przyjąć hipotezę alternatywną. Wynika z tego, że na wartość parametru G ma wpływ zarówno parametr D, jak i rodzaj zastosowanego leku.

22. Porównanie obu modeli

H₀: Brak statystycznie istotnych różnic między zastosowaniem modeli.

H₁: Występują różnice w zastosowanych modelach.

anova(model,model2)

Analysis of Variance Table

Model 1: paramG ~ paramD * lek

Model 2: paramG ~ lek + paramD

Res.Df RSS Df Sum of Sq F Pr(>F)

1 16 109.98

2 17 122.32 -1 -12.337 1.7948 0.1991

Uzyskana w teście p- wartość wynosi 0,199 i jest większa od ustalonego poziomu istotności wynoszącego α= 0,05. Oznacza to, że nie ma podstaw do odrzucenia hipotezy zerowej,
a w związku z tym skonstruowane modele nie różnią się w sposób statystycznie istotny.

Interpretacja analizy kowariancji

Przy wyborze modelu zastosowanego do analizy istnieje zasada, aby wybrany model był jak najmniej złożony („prostszy model jest lepszy”). Postępując zgodnie z tą zasadą wybieramy model 2 i jego wyniki interpretujemy. Według tego testu wcześniej zaproponowana hipoteza zerowa jest błędna i należy przyjąć hipotezę alternatywną. Zgodnie z tą drugą na wartości parametru G mają wpływ zarówno parametr D, jak i rodzaj zastosowanego leku.

(to z lekiem było już widoczne na wykresie w punkcie 4)

Tego nie rozumiem!

## różnica między modelami jest statystycznie nieistotna

## zgodnie z zasadą "prostszy model jest lepszy"

## przyjmujemy model2

## Interpretacja: po uwzględnieniu wpływu parametru D na wartość parametru G

## stwierdzamy statystycznie istotną różnicę między średnimi wartościami parametru G

## w zależności od stosowanego leku

### Uzyskany wniosek jest zatem inny, niż bez uwzględnienia wpływu parametru D

A tego podobno już nie trzeba

20. średnia wartość parametru G, gdy stosowany jest lek B:

w modelu regresji dla pierwszego poziomu zmiennej

## (czyli A) jest przyjmowana wartość zero, a dla drugiego poziomu wartość jeden

18.3600-5.1547*1+0.8275*mean(paramD)

[1] 26.07292

21. średnia wartość parametru G, gdy stosowany jest lek A:

18.3600-5.1547*0+0.8275*mean(paramD)

[1] 31.22762

22. Dorysowanie prostych

• Jako średnia dla leku A

abline(18.3600-5.1547*0,0.8275)

• Jako średnia dla leku B

abline(18.3600-5.1547*1,0.8275,lty=2, col="red")