Labolatorium 3

Karolina Marek

2024-08-06

Test jednorodności wariancji

Dla testów t-Studenta w problemie 2 grup pojawiało się założenie o jednoznacznej wariancji (laboratorium 2).

Dlatego też przedstawimy w tym laboratorium test jednostajnie najmocniejszy wśród testów nieobciążonych oparty na ilorazie wiarygodności, który umożliwi sprawdzenie czy wariancję są jednorodne.

Zacznijmy od przypomnienia definicji rozkładu F.

Rozkład Snedecora-Fishera \(F_{m,n}\)

Niech \(X,Y\) będą niezależnymi zmiennymi losowymi o rozkładzie odpowiednio \(\chi^2_{m}\) oraz \(\chi^2_{n}\). Wtedy \[\frac{\frac{1}{m} X}{\frac{1}{n} Y}\sim F_{m,n}.\] Dodatkowo, \[\frac{Y}{X+Y}\sim Beta\left(\frac{n}{2},\frac{m}{2} \right).\] Opiszemy w tym przypadku test jednostajnie najmocniejszy wśród testów nieobciążonych oparty na ilorazie wiarygodności.

Mamy \(X_1\dots, X_m\sim N(m_x, \sigma^2_X)\) oraz \(Y_1\dots, Y_n\sim N(m_y, \sigma^2_Y)\) - niezależne próby i.i.d.

Testujemy \(H_0: \sigma^2_X=\sigma^2_Y\), przeciwko \(H_1: \sigma^2_X\neq \sigma^2_Y\).

W praktyce, będziemy używać statystyki testowej \[U(X,Y)=\frac{\sum_{i=1}^n\left(Y_i-\bar{Y}_n\right)^2}{\sum_{j=1}^m\left(X_j-\bar{X}_m\right)^2+\sum_{i=1}^n\left(Y_i-\bar{Y}_n\right)^2}\sim Beta\left(\frac{n-1}{2},\frac{m-1}{2} \right).\] Przypomnijmy (szczegóły w wykładzie), że musimy rozważyć dwa przypadki, tj. gdy statystyka \(U(X,Y)\) jest mniejsza od wartości średniej i przypadek odwrotny.

Do praktycznego wyznaczenia wartości krytycznych możemy użyć poniższych warunków: \[F_{B\left(\frac{n-1}{2},\frac{m-1}{2}\right)}\left(D(U)\right)-F_{B\left(\frac{n+1}{2}, \frac{m+1}{2}\right)}\left(D(U)\right)=F_{B\left(\frac{n-1}{2}, \frac{m-1}{2}\right)}\left(U\right)-F_{B\left(\frac{n-1}{2}, \frac{m-1}{2}\right)}\left(U\right),\] \[F_{B\left(\frac{n-1}{2},\frac{m-1}{2}\right)}\left(C(U)\right)-F_{B\left(\frac{n+1}{2}, \frac{m+1}{2}\right)}\left(C(U)\right)=F_{B\left(\frac{n-1}{2}, \frac{m-1}{2}\right)}\left(U\right)-F_{B\left(\frac{n-1}{2}, \frac{m-1}{2}\right)}\left(U\right),\]

Aby wyznaczyć wartości \(D(U)\) oraz \(C(U)\) numerycznie użyjemy metody bisekcji.

B.S.O. załóżmy, że poszukujemy \(D(U)\).

  • W warunku przenosimy wszystko na jedna stronę i przyrównujemy wyrażenie \(\tilde{F}(D)\) do \(0\).

  • Definiujemy \(a_1=0\), \(b_1=(n-1)/(n+m-2)\) oraz \(d_1\) będzie środkiem odcinka, czyli w tym przypadku \(d_1=(a_1+b_1)/2\).

  • Sprawdzamy czy \(\tilde{F}(D)<0\), jeśli tak, to zostawiamy lewy koniec \(a_2=a_1\), a za \(b_2=d_1\).

  • Następnie powtarzamy poprzednie \(2\) kroki określoną ilość razy.

Obliczenie p-wartości: Gdy \(U<\frac{n-1}{n+m-2}\): \[P(X,Y)=F_{B\left(\frac{n-1}{2}, \frac{m-1}{2}\right)}(U)+1-F_{B\left(\frac{n-1}{2}, \frac{m-1}{2}\right)}\left(D(U)\right).\] Gdy \(U>\frac{n-1}{n+m-2}\): \[P(X,Y)=F_{B\left(\frac{n-1}{2}, \frac{m-1}{2}\right)}\left(D(U)\right)+1-F_{B\left(\frac{n-1}{2}, \frac{m-1}{2}\right)}\left(U\right).\]

Zadanie 1 Napisać funkcję testu jednostajnie najmocniejszego wśród testów nieobciążonych dla problemu testowania równości wariancji w modelu Gaussowskim.

calculate_stat_U <- function(X, Y) {
  nx <- length(X)
  ny <- length(Y)
  sum_sq_X <- sum((X - mean(X))^2)
  sum_sq_Y <- sum((Y - mean(Y))^2)
  stat_U <- sum_sq_Y / (sum_sq_Y + sum_sq_X)
  return(list(stat_U = stat_U, nx = nx, ny = ny))
}

calculate_critical_value <- function(nx, ny) {
  return((ny - 1) / (nx + ny - 2))
}

bisection_method <- function(F_tilde, stat_U, critical_value,lower_bound, upper_bound) {
  a <- lower_bound
  b <- upper_bound
  if(stat_U < critical_value) {
    for (i in 1:40) {
    d <- (a + b) / 2
    if (F_tilde(d, stat_U) < 0) {
      b <- d
    } else {
      a <- d
    }
  }
  }
  else{
    for (i in 1:40) {
    d <- (a + b) / 2
    if (F_tilde(d, stat_U) < 0) {
      a <- d
    } else {
      b <- d
    }
  } 
  }
  return((a + b) / 2)
}

calculate_p_value <- function(stat_U, critical_value, d, nx, ny) {
  if (stat_U < critical_value) {
    return(pbeta(stat_U, (ny - 1) / 2, (nx - 1) / 2) + 1 - pbeta(d, (ny - 1) / 2, (nx - 1) / 2))
  } else {
    return(pbeta(d, (ny - 1) / 2, (nx - 1) / 2) + 1 - pbeta(stat_U, (ny - 1) / 2, (nx - 1) / 2))
  }
}

# Główna funkcja
variance_test <- function(X, Y) {
  results <- calculate_stat_U(X, Y)
  stat_U <- results$stat_U
  nx <- results$nx
  ny <- results$ny

  critical_value <- calculate_critical_value(nx, ny)

  F_tilde <- function(x, t) {
    pbeta(x, (ny - 1) / 2, (nx - 1) / 2) - pbeta(x, (ny + 1) / 2, (nx - 1) / 2) - pbeta(t, (ny - 1) / 2, (nx - 1) / 2) + pbeta(t, (ny + 1) / 2, (nx - 1) / 2)
  }

  if (stat_U < critical_value) {
    d <- bisection_method(F_tilde, stat_U, critical_value, critical_value, 1)
  } else {
    d <- bisection_method(F_tilde, stat_U,critical_value, 0, critical_value)
  }

  p_value <- calculate_p_value(stat_U, critical_value, d, nx, ny)
  return(p_value)
}

Zadanie 1 c.d. Przetestować jednorodność wariancji na danych z Baza 2.2.

library(readxl)
Baza<- read_excel("Baza 2.2.xlsx")
X<-na.exclude(Baza$"Zdrowi")
Y<-na.exclude(Baza$"Chorzy")
variance_test(X,Y)
## [1] 0.942379

Nie ma podstaw do odrzucenia hipotezy zerowej.

Możemy również użyć wbudowanego testu do badania jednorodności wariancji - var.test.

var.test(X,Y)
## 
##  F test to compare two variances
## 
## data:  X and Y
## F = 0.97595, num df = 39, denom df = 33, p-value = 0.935
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4966397 1.8824165
## sample estimates:
## ratio of variances 
##          0.9759543

Uwaga: Nie jest to jednak test jednostajnie najmocniejszy wśród nieobciążonych. Oparty jest on na statystyce testowej \(𝑍\), która jest opisana w wykładzie, jednak do obliczenia p-wartości nie używa wartości spełniających założenia twierdzenia o teście jednostajnie najmocniejszym wśród nieobciążonych (po prostu podwaja pole na lewo lub na prawo od statystyki testowej w zależności od jej wartości pod gęstością odpowiedniego rozkładu Snedecora-Fishera, czyli odcina ogony o tym samym prawdopodobieństwie).

Zadanie 2

Narysować wykres empirycznej mocy powyższego testu dla \(20\) obserwacji \(X\) z rozkładu \(N(2,10)\) i \(20\) obserwacji \(Y\) z rozkładu \(N(5,s)\), gdzie \(s\in[1,25]\).

moc<-c()
mc<-1000
s<-c(1:25)
for (k in 1:25){
  m<-0
  for (j in 1:mc){
    X<-rnorm(20,2,10)
    Y<-rnorm(20,5,s[k])
    p<-variance_test(X,Y)
    if(p<0.05){
      m = m+1
      }
    }
  moc[k]<-m/mc
  }
plot(moc)

Zadanie

Sporządzić na jednym rysunku wykres empirycznej mocy dla F-testu (czyli używając gotowej funkcji var.test()) i testu jednostajnie najmocniejszego wśród nieobciążonych (można użyć funkcji z pierwszego zadania) dla \(20\) obserwacji \(X\) z rozkładu \(N(2,10)\) i \(70\) obserwacji \(Y\) z rozkładu \(N(5,s)\), gdzie \(s\in[1,25]\). Czy wykresy różnią się znacząco od siebie?

Przetestuj powyższym kod dla takiej samej liczności prób.

Uwaga: Rozkład Beta jest rozkładem symetrycznym, dlatego też, gdy próby będą różnej wielkości to będzie bardziej widoczna różnica pomiędzy testami.

Test drugi (var.test()) traci moc po prawej stronie - jest obciążony (oszukuje przy dużej dysproporcji pomiędzy próbami).

Zadanie

Pokazać obciążenie gotowej funkcji var.test(). Wygeneruj dwie próby z rozkładu normalnego \(X\sim N(2,5)\) oraz \(Y\sim N(5,s)\), gdzie \(s\in[4,10]\). Aby pokazać obciążenie musi byś dysproporcja w licznościach prób. Narysuj wykres mocy empirycznej i zaznacz poziomą prostą poziom istotności.

Uzyskujemy moce niższe niż założony poziom istotności przy prawdziwości hipotezy zerowej.

Natomiast test var.test jest używany w praktyce, ponieważ jest testem zgodnym, czyli dla \(n\rightarrow\infty\) asymptotycznie będzie zbiegać do testu nieobciążonego.