ZADANIE 4

Kristína Mečiarová


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