Глава 5: Tests on One or Two Mean Vectors

[5.1]

Да се покаже дека:

\[ T^2=n(\bar{\boldsymbol{y}}-\boldsymbol{\mu}_0)'\boldsymbol{S}^{-1}(\bar{\boldsymbol{y}}-\boldsymbol{\mu}_0)= (\bar{\boldsymbol{y}}-\boldsymbol{\mu}_0)'(\frac{\boldsymbol{S}}{n})^{-1}(\bar{\boldsymbol{y}}-\boldsymbol{\mu}_0) \]

\[ (\bar{\boldsymbol{y}}-\boldsymbol{\mu}_0)'(\frac{\boldsymbol{S}}{n})^{-1}(\bar{\boldsymbol{y}}-\boldsymbol{\mu}_0)=(\bar{\boldsymbol{y}}-\boldsymbol{\mu}_0)'(\frac{\boldsymbol{1}}{n})^{-1}\boldsymbol{S}^{-1}(\bar{\boldsymbol{y}}-\boldsymbol{\mu}_0)=n(\bar{\boldsymbol{y}}-\boldsymbol{\mu}_0)'\boldsymbol{S}^{-1}(\bar{\boldsymbol{y}}-\boldsymbol{\mu}_0)=T^2 \]

[5.5]

Да се покаже дека:

\[ s^2_d=\sum_{i=1}^{n}\frac{(d_i-\bar{d})^2}{n-1}=s^2_y+s^2_x-2s_{yx} \]

\[ s^2_d=\sum_{i=1}^{n}\frac{(d_i-\bar{d})^2}{n-1}=\\ \sum_{i=1}^{n}\frac{(y_i-x_i-\frac{1}{n}\sum_{j=1}^{n}(y_j-x_j))^2}{n-1}=\\ \sum_{i=1}^{n}\frac{(y_i-x_i-\bar{y}+\bar{x})^2}{n-1}=\\ \sum_{i=1}^{n}\frac{(y_i-\bar{y})^2+(x_i-\bar{x})^2-2(y_i-\bar{y})(x_i-\bar{x})}{n-1}=\\ s^2_y+s^2_x-2s_{yx} \]

[5.6]

Да се покаже дека:

\[ T^2=n\bar{\boldsymbol{d}}'\boldsymbol{S}_d^{-1}\bar{\boldsymbol{d}}=\bar{\boldsymbol{d}}'(\frac{\boldsymbol{S}_d}{n})^{-1}\bar{\boldsymbol{d}} \]

\[ \bar{\boldsymbol{d}}'(\frac{\boldsymbol{S}_d}{n})^{-1}\bar{\boldsymbol{d}}=\bar{\boldsymbol{d}}'(\frac{1}{n})^{-1}\boldsymbol{S}_d^{-1}\bar{\boldsymbol{d}}=n\bar{\boldsymbol{d}}'\boldsymbol{S}_d^{-1}\bar{\boldsymbol{d}}=T^2 \]

[5.10]

Да се покаже дека:

\[ T^2=(n_1+n_2)(\boldsymbol{C}\bar{\boldsymbol{y}})'(\boldsymbol{C}\boldsymbol{S}_{pl}\boldsymbol{C}')^{-1}\boldsymbol{C}\bar{\boldsymbol{y}} \sim T^2_{p-1,n_1+n_2-2} \]

\(\frac{\boldsymbol{C}\boldsymbol{S}_{pl}\boldsymbol{C}'}{n_1+n_2}\) e коваријансната матрица на \(\boldsymbol{C}\bar{\boldsymbol{y}}\). Кога ќе се земе тоа предвид, изразот комплетно се совпаѓа со карактеристичната форма на \(T^2\) од претходната задача, па може да се заклучи дека \(T^2\) навистина има распределба \(T^2_{p-1,n_1+n_2-2}\).

[5.16]

Податоците се:

data = c(
  1, 189, 245, 137, 163, 1, 181, 305, 184, 209,
  2, 192, 260, 132, 217, 2, 158, 237, 133, 188,
  3, 217, 276, 141, 192, 3, 184, 300, 166, 231,
  4, 221, 299, 142, 213, 4, 171, 273, 162, 213,
  5, 171, 239, 128, 158, 5, 181, 297, 163, 224,
  6, 192, 262, 147, 173, 6, 181, 308, 160, 223,
  7, 213, 278, 136, 201, 7, 177, 301, 166, 221,
  8, 192, 255, 128, 185, 8, 198, 308, 141, 197,
  9, 170, 244, 128, 192, 9, 180, 286, 146, 214,
  10, 201, 276, 146, 186, 10, 177, 299, 171, 192,
  11, 195, 242, 128, 192, 11, 176, 317, 166, 213,
  12, 205, 263, 147, 192, 12, 192, 312, 166, 209,
  13, 180, 252, 121, 167, 13, 176, 285, 141, 200,
  14, 192, 283, 138, 183, 14, 169, 287, 162, 214,
  15, 200, 294, 138, 188, 15, 164, 265, 147, 192,
  16, 192, 277, 150, 177, 16, 181, 308, 157, 204,
  17, 200, 287, 136, 173, 17, 192, 276, 154, 209,
  18, 181, 255, 146, 183, 18, 181, 278, 149, 235,
  19, 192, 287, 141, 198, 19, 175, 271, 140, 192,
  20, 0, 0, 0, 0, 20, 197, 303, 170, 205
)

y1 = data[seq(2, length(data) - 10, 10)]
y2 = data[seq(3, length(data) - 10, 10)]
y3 = data[seq(4, length(data) - 10, 10)]
y4 = data[seq(5, length(data) - 10, 10)]

ole = data.frame(y1, y2, y3, y4)
ole

y1 = data[seq(7, length(data), 10)]
y2 = data[seq(8, length(data), 10)]
y3 = data[seq(9, length(data), 10)]
y4 = data[seq(10, length(data), 10)]

car = data.frame(y1, y2, y3, y4)
car

(a)

Да се тестира хипотезата \(H_0:\boldsymbol{\mu}_1=\boldsymbol{\mu}_2\).

calc_w <- function(df) {
  rows = dim(df)[1]
  cols = dim(df)[2]
  means = colMeans(df)
  w = matrix(data=rep(0, cols * cols), nrow=cols, ncol=cols)
  for (k in 1:rows) {
    row = as.matrix(df[k, ])
    w = w + t(row - means) %*% (row - means)
  }
  return(w)
}

calc_t2_spl <- function(df1, df2) {
  df1_means = colMeans(df1)
  df2_means = colMeans(df2)
  n1 = dim(df1)[1]
  n2 = dim(df2)[1]
  w1 = calc_w(df1)
  w2 = calc_w(df2)
  spl = (w1 + w2) / (n1 + n2 - 2)
  means_diff = df1_means - df2_means
  return(list(
    (t(means_diff) %*% solve(spl) %*% means_diff) * (n1 * n2) / (n1 + n2),
    spl
  ))
}

res = calc_t2_spl(ole, car)
t2 = res[[1]]
spl = res[[2]]
t2

         [,1]
[1,] 133.4873

Оваа \(T^2\) вредност е доволно голема за да ја отфрлиме нултата хипотеза дека средините се исти.

(b)

Да се направи тестот на поединечните променливи.

ole_means = colMeans(ole)
car_means = colMeans(car)
n1 = dim(ole)[1]
n2 = dim(car)[1]

uni_t = (ole_means - car_means) / sqrt((diag(spl) * (n1 + n2) / (n1 * n2)))
uni_t

       y1        y2        y3        y4 
 3.887946 -3.865239 -5.691131 -5.042625 

Вредностите се доволно големи за да се отфрлат нултите хипотези за поединечните променливи т.е. не можеме да заклучиме дека просеците на поединечните променливи се исти.

(c)

Да се пресмета векторот со коефициенти на дискриминантната функција \(\boldsymbol{a}=\boldsymbol{S}_{pl}(\bar{\boldsymbol{y}}_1-\bar{\boldsymbol{y}}_1)\).

a = solve(spl) %*% (ole_means - car_means)
a

         [,1]
y1  0.3452490
y2 -0.1303878
y3 -0.1064338
y4 -0.1433533

(d)

Да се покаже дека \(T^2=t^2(\boldsymbol{a})\).

t_a <- function(a, y1_means, y2_means, n1, n2, spl) {
  return(
    (t(a) %*% y1_means - t(a) %*% y2_means) / 
      sqrt(t(a) %*% spl %*% a * (n1 + n2) / (n1 * n2))
  )
}

t_a(a, ole_means, car_means, n1, n2, spl) ^ 2

         [,1]
[1,] 133.4873

Исто е како и пресметаното под (а).

(e)

Да се најде \(T^2\) користејќи го методот на регресија.

ole["w"] = n2 / (n1 + n2)
car["w"] = n1 / (n1 + n2)
all_data = rbind(ole, car)
all_data

linear_model = lm(w ~ y1 + y2 + y3 + y4, data=all_data)
s = summary(linear_model)
s

Call:
lm(formula = w ~ y1 + y2 + y3 + y4, data = all_data)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.013859 -0.003554  0.000767  0.003522  0.011015 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.231e-01  1.979e-02  26.432  < 2e-16 ***
y1           5.059e-04  8.440e-05   5.994 8.76e-07 ***
y2          -1.911e-04  7.934e-05  -2.408  0.02160 *  
y3          -1.560e-04  1.171e-04  -1.332  0.19163    
y4          -2.101e-04  7.310e-05  -2.874  0.00694 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.006395 on 34 degrees of freedom
Multiple R-squared:  0.783, Adjusted R-squared:  0.7574 
F-statistic: 30.67 on 4 and 34 DF,  p-value: 7.522e-11

r2 = s$r.squared
new_t2 = (n1 + n2 - 2) * (r2 / (1 - r2))
new_t2

[1] 133.4873

Вредноста е иста како и таа пресметана под (а) и под (d).

(f)

Да се пресмета придонесот на секоја променлива во однос на останатите три.

Прикажани се по ред (Ф статистики):

ole = ole[-c(5)]
car = car[-c(5)]

calc_t2_no_ys <- function(df1, df2, omited_ys, t2_total) {
  ni = dim(df1)[1] + dim(df2)[1] - 2
  t_no_ys = calc_t2(df1[-omited_ys], df2[-omited_ys])
  return(((ni - length(df1) + 1) / length(omited_ys)) * (t2_total - t_no_ys) / (ni + t_no_ys))
}

calc_t2_no_ys(ole, car, c(1), t2)

        [,1]
[1,] 35.9336

calc_t2_no_ys(ole, car, c(2), t2)

         [,1]
[1,] 5.799435

calc_t2_no_ys(ole, car, c(3), t2)

         [,1]
[1,] 1.774944

calc_t2_no_ys(ole, car, c(4), t2)

         [,1]
[1,] 8.259241

(g)

Да се пресмета придонесот на последните 2 променливи во однос на првите две.

Прикажана е Ф статистиката:

calc_t2_no_ys(ole, car, c(3, 4), t2)

        [,1]
[1,] 6.08144