V tejto úlohe budeme pracovať s dátovým súborom uscrime
z knižnice SMS data
. Súbor obsahuje 50 pozorovaní a 11
premenných s informáciami o počte výskytov rôznych trestných činov v
rámci USA. Bude nás zaujímať, či sú očakávané počty vrážd, krádeží a
útokov rovnaké v rámci regiónov juh a západ. Pracujeme teda s 2
náhodnými výbermi, z ktorých jeden predstavuje počet vrážd, krádeží a
útokov v južnom regióne a druhý v západnom regióne. Nasledujúca tabuľka
zhŕňa výberové charakteristiky týchto dvoch náhodných výberov.
library(DescTools)
library(kableExtra)
data(uscrime, package = "SMSdata")
data<-uscrime[uscrime$reg == "West" | uscrime$reg == "South",]
data<-transform(data,reg=droplevels(data$reg))
fx <- function(x){
y <- round(c(min(x),max(x), mean(x),median(x), var(x)),2)
names(y) <- c("minimum", "maximum", "priemer", "medián", "rozptyl")
return(y)
}
cX<-rbind(fx(data$murder[data$reg == "South"]),fx(data$robbery[data$reg == "South"]),fx(data$assault[data$reg == "South"]) )
cY<-rbind(fx(data$murder[data$reg == "West"]),fx(data$robbery[data$reg == "West"]),fx(data$assault[data$reg == "West"]) )
tab<-rbind(cX,cY)
rownames(tab)<-c("Vraždy juh", "Krádeže juh", "Útoky juh","Vraždy západ", "Krádeže západ", "Útoky západ")
tab %>%
kbl() %>%
kable_styling()
minimum | maximum | priemer | medián | rozptyl | |
---|---|---|---|---|---|
Vraždy juh | 5.9 | 15.3 | 10.52 | 11.3 | 7.23 |
Krádeže juh | 19.0 | 338.6 | 97.37 | 77.4 | 6089.76 |
Útoky juh | 84.0 | 293.0 | 185.81 | 185.0 | 3431.90 |
Vraždy západ | 3.2 | 12.2 | 6.28 | 5.7 | 8.38 |
Krádeže západ | 20.5 | 206.9 | 89.21 | 71.8 | 3580.32 |
Útoky západ | 43.0 | 226.0 | 135.08 | 132.0 | 3952.08 |
Budeme predpokladať, že máme 2 náhodné výbery z 3-rozmerného normálneho rozdelenia s rôznymi kovariančnými maticami.
layout(matrix(c(0,1,1,0, 2,2,3,3), nrow = 2, byrow = TRUE))
#murder
boxplot(data$murder~data$reg,cex.lab=1.5,cex.main=1.5,cex.axis=1.5,col=c("#CC3333","#990033"),,xlab="",ylab="",names=c("Juh","Západ"))
title(main="Vraždy")
#robbery
boxplot(data$robbery~data$reg,cex.lab=1.5,cex.main=1.5,cex.axis=1.5,col=c("#FFCC66","#CC9933"),xlab="",ylab="",names=c("Juh","Západ"))
title(main="Krádeže")
#assault
boxplot(data$assault~data$reg,cex.lab=1.5,cex.main=1.5,cex.axis=1.5,col=c("#3366CC","#003399"),xlab="",ylab="",names=c("Juh","Západ"))
title(main="Útoky")
Pre náhodné výbery \(\mathbf{X}_1,\dots
\mathbf{X}_{16} \sim
\textit{N}_3(\mathbf{\mu}_X,\mathbb{\Sigma}_X)\) a \(\mathbf{Y}_1,\dots \mathbf{Y}_{13} \sim
\textit{N}_3(\mathbf{\mu}_Y,\mathbb{\Sigma}_Y)\) s troma
odhadovanými parametrami použijeme Hotellingov \(T^2\) test, teda budeme formálne
testovať:
\[ H_0: \mathbf{\mu}_X=\mathbf{\mu}_Y \quad H_A: \mathbf{\mu}_X \neq \mathbf{\mu}_Y .\]
S testovou štatistikou:
\[T^2= \frac{16 \cdot 13 \cdot (16+13−3−1)}{3 \cdot (16+13)^2}(\overline{\mathbf{X}}_{16}−\overline{\mathbf{Y}}_{13})^⊤\mathcal{S}^{−1}(\overline{\mathbf{X}}_{16}−\overline{\mathbf{Y}}_{13}),\]
kde \(\mathcal{S}=\frac{1}{(16+13)}⋅(16
\cdot\mathcal{S}_X + 13 \cdot \mathcal{S}_Y)\) a \(\mathcal{S}_X, \mathcal{S}_Y\) sú výberové
kovariančné matice po rade náhodných výberov \(\mathbf{X}_1,\dots \mathbf{X}_{16}\) a
\(\mathbf{Y}_1,\dots
\mathbf{Y}_{13}\).
Na výstupe pozorujeme p-hodnotu rovnú \(p=0,006335\), čiže zamietame nulovú hypotézu. Vyvrátili sme, že by bol súčasne počet vrážd, krádeží a útokov na juhu a západe USA zhodný.
X<-cbind(data$murder[data$reg == "South"],data$robbery[data$reg == "South"],data$assault[data$reg == "South"])
Y<-cbind(data$murder[data$reg == "West"],data$robbery[data$reg == "West"],data$assault[data$reg == "West"])
(HT<-HotellingsT2Test(X,Y,mu=NULL,var.equal=FALSE))
##
## Hotelling's two sample T2-test
##
## data: X and Y
## T.2 = 5.1886, df1 = 3, df2 = 25, p-value = 0.006335
## alternative hypothesis: true location difference is not equal to c(0,0,0)
V nasledujúcej tabuľke máme intervaly spoľahlivosti bodové odhady rozdielov stredných hodnôt \(\mathbf{\mu}_X\) a \(\mathbf{\mu}_Y\). Pozorujeme, že ani jeden interval spoľahlivosti nepokrýva hodnotu 0, teda aj z duality štatistických testov a intervalových odhadov zamietame nulovú hypotézu.
S<-1/29*(1/16*cov(X)+1/13*cov(Y))
Fquant <- (3/(16+13 -4)) * qf(1 - 0.05, 3, 25)
meanVec <- c(tab[1,3]-tab[4,3],tab[2,3]-tab[5,3],tab[3,3]-tab[6,3])
LB <- diag(3) %*% meanVec - sqrt(Fquant * diag(diag(3) %*% S %*% diag(3)))
UB <- diag(3) %*% meanVec + sqrt(Fquant * diag(diag(3) %*% S %*% diag(3)))
confInt <- data.frame(LB, UB, meanVec, row.names = c("Vraždy","Krádeže","Útoky"))
names(confInt) <- c("Lower Boundary", "Upper Boundary", "Mean Estimate")
confInt %>%
kbl() %>%
kable_styling()
Lower Boundary | Upper Boundary | Mean Estimate | |
---|---|---|---|
Vraždy | 4.123497 | 4.356503 | 4.24 |
Krádeže | 5.310455 | 11.009545 | 8.16 |
Útoky | 48.196672 | 53.263328 | 50.73 |