Test Welscha
W przypadku, gdy wariancja nie jest jednorodna nie możemy używać testu t-Studenta.
Załóżmy, że mamy następująca sytuację: \(X\sim N(m_x, \sigma^2_x)\), \(Y\sim N(m_y, \sigma^2_y)\) i chcemy przetestować \(H_0: m_x=m_y\) przeciwko \(H_1: m_x\neq m_y\).
Wtedy statystyka \[W(X,Y)=\frac{\bar{X}-\bar{Y}}{\sqrt{\frac{1}{m(m-1)}\sum_{i=1}^m(X_i-\bar{X})^2+\frac{1}{n(n-1)}\sum_{i=1}^n(Y_i-\bar{Y})^2}}\] ma rozkład asymptotycznie \(N(0,1)\) dla \(m,n\rightarrow \infty\), przy założeniu \(H_0\).
Aby utrzymać założony poziom istotności będziemy dla mniejszych \(m,n\), użyjemy testu Welscha.
Gdy \(m_x=m_y\), to dla dowolnego \(\epsilon>0\) zachodzi \[\lim_{\min\{m,n\}\rightarrow\infty}P\left(\sup_{x\in \mathbb{R}}\left|F_W(x)-F_{t_{R(x,Y)}}\right|<\epsilon\right)=1,\]
gdzie \(F_W(x)\) jest nieznaną dystrybuantą zmiennej losowej \(W\), a \(F_{t_{R(X,Y)}}\) jest dystrybuantą rozkładu t-Studenta o \(R(X,Y)\) stopniach swobody, gdzie
\[R(X,Y)=\frac{\left(\frac{m}{m-1}\sum_{i=1}^m(X_i-\bar{X})^2+\frac{m}{n-1}\sum_{i=1}^n(Y_i-\bar{Y})^2\right)(m-1)(n-1)}{\left(\frac{n}{m-1}\sum_{i=1}^m(X_i-\bar{X})^2\right)(n-1)+\left(\frac{m}{n-1}\sum_{i=1}^n(Y_i-\bar{Y})^2\right)(m-1)}.\]
Zadanie 1
W grupie 20 sportowców oraz w grupie kontrolnej 20 osób wykonano pomiar tętna spoczynkowego. Sprawdzić założenie jednorodności wariancji. Przetestować hipotezę o tym, że tętno spoczynkowe u sportowców jest niższe niż u osób nie uprawiających zawodowo sportu (2 sposoby: oprogramować i gotowiec). Dane są w pliku Baza 4.1.
Przeprowadzamy test \(H_0: m_x<m_y\), przeciwko \(H_1: m_x\geq m_y\).
Wtedy, \(\phi(X,Y)=1\) dla \(W(X,Y)>c_{\alpha}(X,Y)\), gdzie \(c_{\alpha}=F^{-1}_{t_{R(X,Y)}}(1-\alpha)\).
P-wartość: \(p(X,Y)=1-F_{t_{R(X,Y)}}\left(W(X,Y)\right)\)
library(readxl)
Baza<- read_excel("Baza 4.1.xlsx")
X<-Baza$"Sportowcy"
Y<-Baza$"Kontrolna"
m<-length(X)
n<-length(Y)
var.test(X,Y) # test na jednorodność wariancji ##
## F test to compare two variances
##
## data: X and Y
## F = 0.37391, num df = 19, denom df = 19, p-value = 0.03789
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.1479998 0.9446760
## sample estimates:
## ratio of variances
## 0.3739143##
## Welch Two Sample t-test
##
## data: X and Y
## t = -2.7051, df = 31.466, p-value = 0.005464
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
## -Inf -3.249356
## sample estimates:
## mean of x mean of y
## 60.25 68.95W_stat<-(mean(Y)-mean(X))/sqrt(sum((X-mean(X))^2)/(m*(m-1))+sum((Y-mean(Y))^2)/(n*(n-1)) )
R<-((n-1)*(m-1)*(sum((X-mean(X))^2)*(n/(m-1))+sum((Y-mean(Y))^2)*(m/(n-1)))^2)/((n-1)*(sum((X-mean(X))^2)*(n/(m-1)))^2+(m-1)*(sum((Y-mean(Y))^2)*(m/(n-1)))^2)
p<-1-pt(W_stat,R)
p## [1] 0.005463609Odrzucamy hipotezę zerową na rzecz alternatywy.
Zadanie
W grupie \(20\) sportowców oraz w grupie kontrolnej \(20\) osób wykonano pomiar tętna spoczynkowego. Przetestować hipotezę o tym, że tętno spoczynkowe u sportowców jest takie samo jak u osób nie uprawiających zawodowo sportu. Dane są w pliku Baza 4.1. Porównać do testu t-Studenta wykonywanego przy założeniu jednorodności wariancji.
Test niezależności dla rozkładów normalnych
W teście t-Studenta zakładamy nie tylko jednorodności wariancji, ale również niezależność obserwacji.
Jeżeli założymy, że obserwujemy pary \((X_1, Y_1), \dots, (X_n, Y_n)\), to istnieje test jednostajnie najmocniejszy wśród testów nieobciążonych, który bada czy obserwację te są niezależne.
Dla rozkładów normalnych niezależność jest równoważna zerowaniu się współczynnika korelacji.
Testujemy zatem niezależność \(H_0: \rho=0\) przeciwko \(H_1: \rho\neq0\).
Wobec tego, niech statystyka testowa (wynikająca z test opartego na ilorazie wiarygodności) będzie następującej postaci: \[U(X,Y)=\hat{\rho}(X,Y)=\frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})}{\sqrt{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})}},\] ma rozkład symetryczny wokół zera.
Zatem, test \(\phi(X,Y)=1\) dla \(|\hat{\rho}(X,Y)|>c_{\alpha}\), gdzie \(c_{\alpha}=2F^{-1}_{B\left(\frac{n}{2}-1,\frac{n}{2}-1\right)}\left(1-\frac{\alpha}{2}\right)\) jest testem jednostajnie najmocniejszym wśród testów nieobciążonych.
P-wartość: \[p(X,Y)=2\left(1-F_{B\left(\frac{n}{2}-1, \frac{n}{2}-1\right)}\left(\frac{1+|\hat{\rho}(X,Y)|}{2}\right)\right).\]
Zadanie 2
\(20\) pacjentom zmierzono ciśnienie tętnicze krwi (skurczowe i rozkurczowe). Na podstawie zebranych danych przetestować hipotezę o niezależności ciśnienia skurczowego od rozkurczowego. Dane są w pliku Baza 4.2.
library(readxl)
Baza<- read_excel("Baza 4.2.xlsx")
X<-Baza$"Skurczowe"
Y<-Baza$"Rozkurczowe"
n<-length(X)
plot(X,Y)
##
## Pearson's product-moment correlation
##
## data: X and Y
## t = 6.7271, df = 18, p-value = 2.633e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6445724 0.9374619
## sample estimates:
## cor
## 0.8458311rho<-sum((X-mean(X))*(Y-mean(Y)))/sqrt(sum((X-mean(X))^2)*sum((Y-mean(Y))^2))
p<-2*(1-pbeta((1+abs(rho))/2,n/2-1,n/2-1))
rho## [1] 0.8458311## [1] 2.632736e-06Zatem odrzucamy hipotezę zerową na rzecz alternatywy.
Możemy również użyć gotowej funkcji cor.test().
##
## Pearson's product-moment correlation
##
## data: X and Y
## t = 6.7271, df = 18, p-value = 2.633e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6445724 0.9374619
## sample estimates:
## cor
## 0.8458311Zadanie 3
Rozważmy 20 elementową próbę prostą wektora \((X,Y)\), gdzie \(Y=aX+Z\),\(X\sim N(10,5)\), \(Z\sim N(0,s)\) oraz niech X, Z będą niezależne. Sporządzić wykres mocy empirycznej testu niezależności \((X,Y)\) w zależności od \((a,s)\), gdzie \(a\in[-2,2]\) oraz \(s\in[1,20]\).
mc<- 1000
n <- 20
alpha<-0.05
Moc<- matrix(0,nrow=20,ncol=20)
r <- seq(-2,2, length.out=20)
for (j in 1:20){
for (k in 1:20){
a<- r[j]
s<-k
M <- 0
for (i in 1:mc){
X <- rnorm(n,10,5)
Y<- a*X+rnorm(n,0,s)
rho<-sum((X-mean(X))*(Y-mean(Y)))/sqrt(sum((X-mean(X))^2)*sum((Y-mean(Y))^2))
p=2*(1-pbeta((1+abs(rho))/2,n/2-1,n/2-1))
if(p< alpha){
M=M+1
}
}
Moc[k,j]<-M/mc
}
}
heatmap(Moc, Rowv = NA, Colv = NA, scale = "none",
labRow = 1:20, labCol = round(r, 2), xlab = "a", ylab = "s",
main = "Moc empiryczna testu niezależności w zależności od a i s")
Zmienne \(X\), \(Y\) będą niezależne, gdy \(a=0\).
ANOVA
Rozważmy przypadek, gdy chcemy sprawdzić równość średnich w więcej niż 2 grupach. Niech \[X_{1,1}, \dots, X_{1, n_1}\sim N(m_1, \sigma^2)\\ \vdots\\ X_{k,1}, \dots, X_{k, n_k}\sim N(m_k, \sigma^2)\] będą niezależnymi próbami prostymi.
Testujemy \(H_0: m_1=\dots=m_k\) przeciwko \(H_1: \exists i,j\in\{1, \dots, k\}\quad m_j\neq m_i\). Oznaczmy \(n=n_1, \dots, n_k=n\).
Test jest oparty na ilorazie wiarygodności.
Otrzymujemy statystykę testową postaci: \[F(\mathbf{X})=\frac{\frac{1}{k-1}\sum_{i=1}^k n_i\left(\bar{X}_{i, n_i}-\bar{X}_n\right)^2}{\frac{1}{n-k}\sum_{i=1}^k \sum_{j=1}^{n_i}\left(X_{i,j}-\bar{X}_{i, n_i}\right)^2},\] która przy prawdziwości hipotezy zerowej ma rozkład Snedecora-Fishera \(F_{k-1, n-k}\).
Otrzymujemy test postaci: \(\phi(\mathbf{X})=1\) dla \(F(\mathbf{X})>c_{\alpha}\), gdzie \(c_{\alpha}=F_{k-1, n-k}^{-1}(1-\alpha)\).
P-wartość: \(p(\mathbf{X})=1-F_{k-1, n-k}\left(F\mathbf(X)\right)\).
Zadanie 4
W 4 klasach przeprowadzono ten sam test z matematyki. Wyniki procentowe zapisane są w pliku Baza 4.3. Przetestować hipotezę o braku różnicy między klasami.
library(readxl)
Baza<- read_excel("Baza 4.3.xlsx")
X<-Baza$"Wynik"
Y<-Baza$"Klasa"
ANOV_func <- function(X, Y) {
n <- length(X)
unique_groups <- unique(Y)
k <- length(unique_groups)
if (k < 2) {
stop("Musi być co najmniej 2 różne grupy.")
}
N <- numeric(k)
Z <- numeric(k)
W <- numeric(n)
for (i in 1:k) {
group <- unique_groups[i]
group_indices <- which(Y == group)
N[i] <- length(group_indices)
Z[i] <- mean(X[group_indices])
W[group_indices] <- X[group_indices] - Z[i]
}
SSB <- sum(N * (Z - mean(X))^2)
SSW <- sum(W^2)
F_stat <- (SSB / (k - 1)) / (SSW / (n - k))
p <- 1 - pf(F_stat, k - 1, n - k)
return(p)
}
ANOV_func(X,Y)## [1] 0.3322455Możemy również użyć gotowej funkcji aov i wyświetlić
wyniki za pomocą summary(aov()).
Baza2<- read_excel("Baza 4.3.xlsx")
Baza2$Klasa <- as.factor(Baza2$Klasa)
Test<-aov(Wynik ~ Klasa, data = Baza2)
summary(Test)## Df Sum Sq Mean Sq F value Pr(>F)
## Klasa 3 1127 375.8 1.151 0.332
## Residuals 102 33296 326.4Nie ma podstaw do odrzucenia hipotezy zerowej.