Labolatorium 5

Karolina Marek

2024-08-06

Testy sympotowtycznie zgodne

Przechodzimy do klasy testów które są asymptotycznie zgodnie.

W przypadku, gdy nie potrafimy znaleźć rozkładu statystyki testowej, spróbujemy wyznaczyć asymptotyczny rozkład i skonstruować test asymptotycznie zgodny.

gdzie \(\labda(\mathbf{X})\) - iloraz największej wiarygodności.

G-test

Jest to test zgodności dla rozkładu wielowymiarowego.

Niech \(X_1, \dots, X_n\) będzie próbą prostą z rozkładu dyskretnego na zbiorze \(\{a_1, \dots, a_k\}\) oraz \[p_j=P(X_i=a_j)>0, \quad j=1,\dots, k\] będą nieznanymi parametrami. Zatem \(\Theta=\{\theta\in(p_1, \dots, p_{k-1}):\quad p_j>0,\quad \sum_{j=1}^{k-1}p_j<1\}\)

Testujemy \(H_0: p_1=\pi_1, \dots, p_k=\pi_k\) przeciwko \(H_1: \exists j\in\{1, \dots, k\} p_j\neq \pi_j\).

Iloraz wiarygodności jest następującej postaci

\[\lambda(N)=\left(\frac{N_1}{n\pi_1}\right)^{N_1}\cdot\dots\cdot \left(\frac{N_k}{n\pi_k}\right)^{N_k},\] gdzie \(N_j(X)=\sum_{i=1}^n\mathbb{I}_{\{a_j\}}(X_i)\).

Korzystając z twierdzenia, otrzymujemy \[G(N)=\sum_{j=1}^k 2 N_j \log\left(\frac{N_j}{n\pi_j}\right)\sim \chi^2_{k-1}.\]

Wobec tego test ma następującą postać \(\phi(N)=1\), gdy \(G(N)>c_{\alpha}\), gdzie \(c_{\alpha}=F^{-1}_{\chi^2_{k-1}}(1-\alpha)\)

P-wartość: \(p(N)=1-F_{\chi^2_{k-1}}\left(G(N)\right)\).

Zadanie 1

Rzucono 40 razy kostką do gry (wyniki w Baza 5.1). Na podstawie uzyskanych wyników przetestować czy kostka jest rzetelna.

library(readxl)
library(DescTools)
Baza<- read_excel("Baza 5.1.xlsx")
X<-Baza$"Wyniki"
n<-length(X)

N <- table(factor(X, levels = 1:6))

P <- rep(1/6, 6)

G <- sum(ifelse(N != 0, 2 * N * log(N / (n * P)), 0))

p_value <- 1 - pchisq(G, df = 5)
p_value
## [1] 0.4462871

Obliczenie p-wartości wykorzystując Monte Carlo

Asymptotyka testu chi kwadrat czasem jest słaba dla małej liczby obserwacji, co powoduje że poziom istotności może nie być zachowany lub moc testu może spaść. W takim przypadku można zastosować metodę Monte Carlo do obliczenia dokładnej p-wartości.

  1. Generujemy mc prób prostych \(X_1, \dots, X_n\) z rozkładu \(H_0\), tj. \(P(X_i=a_j)=\pi_j\).

  2. Dla każdej próby z kroku 1 obliczamy wartość statystyki \(G(N)\). Otrzymujemy próbę prostą \(G_1, \dots, G_{mc}\).

  3. W calu obliczenia p-wartości, obliczamy statystykę \(G(N)\) z oryginalnej (niewygenerowanej) próby. Wtedy \[\hat{p}=\frac{1}{mc}\sum_{i=1}^{mc}\mathbb{I}\left(G_i\geq G(N)\right).\]

Zadanie 1 - c.d.

set.seed(123)
MC <- 100000
GY <- numeric(MC)


for (k in 1:MC) {
  Y <- sample(1:6, n, replace = TRUE)
  NY <- table(factor(Y, levels = 1:6))
  GY[k] <- sum(ifelse(NY != 0, 2 * NY * log(NY / (n * P)), 0))
}

pMC <- mean(GY >= G)
pMC
## [1] 0.4708

Klasyczny test Chi-kwadrat

W praktyce często używa się testu zgodności \(\chi^2\). Nie różni się on zbytnio od G-testu.

Testujemy \(H_0: p_1=\pi_1, \dots, p_k=\pi_k\), przeciwko \(H_1: \exists j\in\{1, \dots, k\} p_j\neq \pi_j\).

Niech \[Q(N)=\sum_{j=1}^k\frac{(N_j-n\pi_j)^2}{2n\pi_j}.\]

Można pokazać, że przy prawdziwości hipotezy zerowej \(\frac{N_j}{n\pi_j}\xrightarrow{n\rightarrow\infty}1\), statystyki \(G(N)\) i \(Q(N)\) mają takie same rozkłady asymptotyczne \(\sim\chi^2_{k-1}\).

Mamy zatem test: \(\phi(N)=1\) dla \(Q(N)>c_{\alpha}\), gdzie \(c_{\alpha}=F^{-1}_{\chi^2_{k-1}}(1-\alpha)\).

P-wartość: \(p(N)=1-F_{\chi^2_{k-1}}\left(Q(N)\right)\).

Interpretacja testu \(\chi^2\)

  • \(\frac{(N_j-n\pi_j)}{n\pi_j}\) - kwadrat przeskalowanej odległości liczności obserwowanej od oczekiwanej (przy prawdziwości hipotezy zerowej).

Zadanie 2

Przeprowadź test \(\chi^2\) rzetelności kostki dla danych Baza 5.1.

library(readxl)
Baza<- read_excel("Baza 5.1.xlsx")
X<-Baza$"Wyniki"
n<-length(X)

N <- table(factor(X, levels = 1:6))
P <- rep(1/6, 6)
Q= sum((N-n*P)^2/(n*P))
p<-1-pchisq(Q,5)
p
## [1] 0.4535778

Możemy również użyć gotowej funkcji chisq.test()

chisq.test(N,p=P)
## 
##  Chi-squared test for given probabilities
## 
## data:  N
## X-squared = 4.7, df = 5, p-value = 0.4536

Za pomocą funkcji chisq.test możemy również obliczyć dokładną (nie asymptotyczną) p-wartość metodą Monte Carlo.

chisq.test(N,p=P,simulate.p.value = TRUE,B=100000)
## 
##  Chi-squared test for given probabilities with simulated p-value (based
##  on 1e+05 replicates)
## 
## data:  N
## X-squared = 4.7, df = NA, p-value = 0.468

Zadanie - obliczenie p-wartości korzystając z metedy Monte Carlo

Oblicz dokładną (nie asymptotyczną) p-wartość testu \(\chi^2\) metodą Monte Carlo nie używając funkcji chisq.test.

Test niezależności \(\chi^2\)

Rozpatrzymy teraz przypadek złożonej hipotezy dla testu zgodności w rozkładzie wielomianowym.

Dla takiej hipotezy będziemy zakładać, że nieznane prawdopodobieństwa komórek dadzą się zapisac jako pewne nieznane funkcje \(d\) zmiennych gdzie \(d<k-1\).

Niech \((X_1, Y_1), \dots, (X_n, Y_n)\) będzie próbą prostą z rozkładu dyskretnego na zbiorze \(\{a_1, \dots, a_k\} \times \{b_1, \dots, b_m\}\). Niech \(p_{ij}=P\left(X_1= a_i, Y_1=b_j\right)\) dla \(i=1, \dots, k\), \(j=1, \dots, m\).

Wtedy

\(p_{i\cdot}=(p_{1,\cdot}, \dots, p_{k\cdot})\) - rozkład brzegowy X,

\(p_{\cdot j}=(p_{\cdot 1}, \dots, p_{\cdot m})\) - rozkład brzegowy Y.

Testujemy: \[H_0: \forall (i,j)\in\{1, \dots, k\}\times \{1, \dots, m\} \quad p_{ij}=p_{i\cdot}p_{\cdot j}\]

\[H_1: \exists (i,j)\in\{1, \dots, k\}\times \{1, \dots, m\} \quad p_{ij}\neq p_{i\cdot}p_{\cdot j}\]

Niech \(\mathbf{N}={N_{11}, \dots, N_{km}}\) - wektor liczności, mający rozkład wielomianowy \(W(n, p_{11}, \dots, p_{km})\).

Sformułujmy hipotezę tak, aby móc skorzystać z twierdzenia. \[H_0: \exists \mathbf{\bar{p}}=(p_{1\cdot}, \dots, p_{k-1\cdot}, p_{\cdot 1}, \cdot, p_{\cdot, m-1})\in\mathbb{R}^{(k-1)+(m-1)}\quad\\ (p_{11}, \dots, p_{km})=\left(g_{11}(p_{1\cdot}, \dots, p_{k-1\cdot}, p_{\cdot 1}, \cdot, p_{\cdot, m-1}), \dots, g_{km}(p_{1\cdot}, \dots, p_{k-1\cdot}, p_{\cdot 1}, \cdot, p_{\cdot, m-1})\right)\] \[H_1: \forall \mathbf{\bar{p}}=(p_{1\cdot}, \dots, p_{k-1\cdot}, p_{\cdot 1}, \cdot, p_{\cdot, m-1})\in\mathbb{R}^{(k-1)+(m-1)}\quad\\ (p_{11}, \dots, p_{km})\neq\left(g_{11}(p_{1\cdot}, \dots, p_{k-1\cdot}, p_{\cdot 1}, \cdot, p_{\cdot, m-1}), \dots, g_{km}(p_{1\cdot}, \dots, p_{k-1\cdot}, p_{\cdot 1}, \cdot, p_{\cdot, m-1})\right)\]

Estymatorem największej wiarygodności wektora \(\mathbf{\bar{p}}\) jest frakcja \[\mathbf{\hat{p}}=\left(\frac{1}{n}\sum_{j=1}^mN_{1j}, \dots, \frac{1}{n}\sum_{j=1}^mN_{kj}, \frac{1}{n}\sum_{j=1}^kN_{i1}, \dots, \frac{1}{n}\sum_{j=1}^kN_{im}\right).\]

Statystyka testowa: \[G(\mathbf{N})=\sum_{i=1}^k\sum_{j=1}^m\frac{\left(N_{ij}-\frac{1}{n}\left(\sum_{i=1}^kN_{ij}\right)\left(\sum_{j=1}^mN_{ij}\right)\right)^2}{\frac{1}{n}\left(\sum_{i=1}^kN_{ij}\right)\left(\sum_{j=1}^mN_{ij}\right)}\] Wtedy \(\phi(\mathbf{N})=1\) dla \(Q(\mathbf{N})>c_{\alpha}\), gdzie \(c_{\alpha}=F^{-1}_{\chi^2_{km-(k+m-2)-1}}(1-\alpha)\).

P-wartość: \(p(\mathbf{N})=1-F_{\chi^2_{km-(k+m-2)-1}}\left(Q(\mathbf{N})\right)\).

Zadanie 3

W latach sześćdziesiątych przeprowadzono badania wśród amerykańskiej młodzieży szkolnej mające wyjaśnić, czy istnieje związek pomiędzy poglądami politycznymi a paleniem marihuany (wyniki w Baza 5.2). Zweryfikować hipotezę o niezależności stosunku do palenia marihuany (1-nigdy, 2-okazjonalnie, 3-często) od poglądów politycznych (1-liberalne, 2-konserwatywne, 3- inne).

library(readxl)
Baza<- read_excel("Baza 5.2.xlsx")
X<-Baza$"Marichuana"
Y<-Baza$"Poglądy"
table(X,Y)
##    Y
## X     1   2   3
##   1 479 214 172
##   2 173  47  45
##   3 119  15  85
n<-length(X)
N<-table(X,Y)
m<-nrow(N)
k<-ncol(N)
Q<-0
for (i in 1:m){
  for (j in 1:k){
    Q=Q+(N[i,j]-sum(N[i,])*sum(N[,j])/n)^2/(sum(N[i,])*sum(N[,j])/n)
    }
  }
p<-1-pchisq(Q,k*m-(k+m-2)-1)
p
## [1] 3.043121e-13

Możemy również użyć gotowej funkcji chisq.test `

chisq.test(table(X, Y))
## 
##  Pearson's Chi-squared test
## 
## data:  table(X, Y)
## X-squared = 64.654, df = 4, p-value = 3.043e-13

Test zgodności Chi-kwadrat hipoteza złożona

Zadanie 4

Cechy dziedziczne są przekazywane za pomocą genów występujących parami. W najprostszym przypadku, każdy gen może mieć dwie formy (allele), które będziemy oznaczać odpowiednio literami A i a. Oznaczmy przez \(P(A)=\theta\) nieznane prawdopodobieństwo występowania genu A w populacji. Z teorii Mendla losowego tworzenia się genotypów potomstwa otrzymujemy następujący rozkład prawdopodobieństwa genotypów populacji w stanie równowagi (równowaga Hardy’ego Weinberga)

Zaobserwowano następujące liczności genotypów:

Czy na poziomie 0.05 są podstawy do odrzucenia hipotezy równowagi H-W?

N<-c(110,235,155)
n<-sum(N)
theta<-(2*N[1]+N[2])/(2*n)
P<-c(theta^2,2*theta*(1-theta),(1-theta)^2)
Q<-0
for (j in 1:3){
  Q=Q+(N[j]- (n*P[j]))^2/(n*P[j])
  }
p<-1-pchisq(Q,1)
p
## [1] 0.2420025

Zadanie

Wykażmy, że estymator największej wiarygodności \(\theta\) jest \(\hat{\theta}=\frac{2n_{AA}+n_{Aa}}{2n}\), gdzie \(n=n_{AA}+n_{Aa}+n_{aa}\).