Практическая работа №3

Проверка гипотез однородности

Глушков Егор Александрович, гр. 20.М04-мм

---

Данные (addicts.xls). Варианты метрической переменной (variable), категориальной с двумя градациями (factor.2), категориальной с четырьмя градациями (factor.4) представлены в Таблице 2.

• Проверить гипотезу о равенстве дисперсий двух выборок и в соответствии с выводом применить критерий Стьюдента для проверки равенства средних. Использовать вариант группирующей переменной factor.2 (Табл. 2).

• Применить однофакторный дисперсионный анализ в случае фактора с четырьмя градациями и множественные сравнения с разными поправками. Проверить гипотезу о равенстве дисперсий.

• Повторить обработки с применением непараметрических аналогов.

• Для первых двух зависимых переменных (Табл. 3, данные dataNF.xls) проверить однородность изменений во времени по критерию Стьюдента для зависимых выборок и по ранговому критерию Вилкоксона.

• Для зависимых переменных (Табл. 3, данные dataNF.xls) с факторами “PRCOD.1” и “SEX.1” выполнить ANOVA Repeated Measures. Проверить значимость факторов “PRCOD.1” и “SEX.1” времени и эффектов взаимодействия.

Анализ независимых выборок

Переменные (вариант 12):

• bdi – variable – оценка депрессии
• se – factor.2 – использование успокоительных
• educat – factor.4 – образование

Исследуем переменные на наличие пропусков. Выделим нужные переменные

data <- na.omit(addicts[ , c("bdi", "se", "educat")])
summary(data)

      bdi              se             educat     
 Min.   : 2.00   Min.   :0.0000   Min.   :1.000  
 1st Qu.:15.00   1st Qu.:0.0000   1st Qu.:2.000  
 Median :20.00   Median :0.0000   Median :2.000  
 Mean   :21.02   Mean   :0.2355   Mean   :2.087  
 3rd Qu.:27.00   3rd Qu.:0.0000   3rd Qu.:2.000  
 Max.   :48.00   Max.   :1.0000   Max.   :4.000  

summary(as.factor(data$se))

  0   1 
211  65 

summary(as.factor(data$educat))

  1   2   3   4 
 21 219  27   9 

library('lawstat')

DescriptiveStat <- function(X, group)
{
  mm. <- tapply(X, group, function(x) mean(x, na.rm=TRUE)); mm.
  sd. <- tapply(X, group, function(x) sd(x, na.rm=TRUE)); sd.
  nn. <- tapply(X, group, function(x) length(na.omit(x))); nn.
  err. <- sd./sqrt(nn.); err.
  list(mm=mm., sd=sd., nn=nn., err=err.)
}

Fig <- function(x)
{
  hist(x, freq=FALSE)
  f1 <- function(x) dnorm(x, mean(x, na.rm=TRUE), sd(x, na.rm=TRUE))
  curve(f1, min(x), max(x), col=2, add=TRUE)
  title(sub=paste("p.Shapiro", format(shapiro.test(x)$p.value, 4, 2), sep="="))
}

Sentence <- function(mm, err, nn, p.T)
{
  A1 <- paste(paste(format(mm, digits=3, nsmall=2), format(err, digits=2, nsmall=2), sep="±"), nn, sep="/")
  A1. <- paste("The means of two groups are", paste(A1, collapse=", "))
  A2. <- ifelse(p.T>0.05, "difference is insignificant", "difference is significant" )
  A3. <- paste("p", format(p.T,3,3),sep="=")
  paste(c(A1., A2., A3.), collapse=", ")
}

Фактор с двумя градациями

Переменная se – группирующая, имеет 2 градации, в роли метрической переменной – bdi

table(data$se) # 2 градации, группирующая переменная

  0   1 
211  65 

df <- data.frame(group=data$se, X=data$bdi)

p.Sh <- with(df, tapply(X, group, function(x)shapiro.test(x)$p.value)); p.Sh

        0         1 
0.5125098 0.4573447 

Согласие с нормальным распределением для обеих групп (градаций) не отвергается в соответствии с критерием Шапиро-Уилка с \(p.value\) 0.51 и 0.457 для обеих групп.

boxplot(X~group, df)

p.F <- var.test(X~group, df)$p.value; p.F

[1] 0.7522756

В соответствии с критерием Фишера гипотеза о равенстве дисперсий в двух групп не отвергается (\(p.value = 0.75\)), значит, можно использовать критерий Стьюдента (с параметром о равенстве дисперсий).

p.T <- t.test(X~group, df, var.equal=TRUE)$p.value; p.T

[1] 0.2282577

По критерию Стьюдента гипотеза о равенстве средних также не отвергается (\(p.value = 0.228\)).

op <- par(mfrow=c(1,2))
Fig(df$X[df$group==0])
Fig(df$X[df$group==1])

par(op)
df.2 <- df

L <- DescriptiveStat(df$X, df$group)
Sentence(L$mm, L$err, L$nn, p.T)

[1] "The means of two groups are 21.37±0.60/211, 19.88±1.11/65, difference is insignificant, p=0.228"

Фактор с четырьмя градациями

Переменная educat – группирующая, имеет 4 градации, в роли метрической переменной – bdi

df <- data.frame(group=as.factor(data$educat), X=as.numeric(data$bdi))
table(df$group)

  1   2   3   4 
 21 219  27   9 

name.gr <- "educat"
name.x <- "bdi"

bartlett.test(X~group, df)

    Bartlett test of homogeneity of variances

data:  X by group
Bartlett's K-squared = 3.7037, df = 3, p-value = 0.2953

Критерий Барлетта показывает, что гипотеза о равенстве дисперсий всех выборок не отвергается с \(p.value = 0.2953\).

with(df,levene.test(X, group))

    Modified robust Brown-Forsythe Levene-type test based on the absolute deviations from
    the median

data:  X
Test Statistic = 1.4306, p-value = 0.2341

Аналогично критерий Левена показывает, что гипотеза о равенстве дисперсий всех выборок не отвергается с p_value = 0.2341.

---

Однофакторный дисперсионный анализ

ao <- aov(X~group, df)
summary(ao)

             Df Sum Sq Mean Sq F value Pr(>F)
group         3    457  152.31   2.012  0.113
Residuals   272  20593   75.71               

boxplot(X~group, xlab=name.gr, ylab=name.x, data=df)

L <- DescriptiveStat(df$X, df$group); L

$mm
       1        2        3        4 
22.00000 21.37900 17.14815 21.66667 

$sd
        1         2         3         4 
 9.137833  8.497583  8.401940 12.971122 

$nn
  1   2   3   4 
 21 219  27   9 

$err
        1         2         3         4 
1.9940387 0.5742134 1.6169541 4.3237073 

Заметим, что средние в каждой из групп примерно равны, критерии выше лишь не опровергают это. При этом различия могут быть объяснены случайностью или малым числом наблюдений в некоторых группах (о чем свидетельствуют параметры nn, sd, err, например, для группы “4”).

---

Множественные сравнения

library(agricolae)
library(multcomp)

ao <- aov(X~group, df)
out <- LSD.test(ao,"group", p.adj="none", group=FALSE); out

$statistics
   MSerror  Df     Mean       CV
  75.70938 272 21.02174 41.39103

$parameters
        test p.ajusted name.t ntr alpha
  Fisher-LSD      none  group   4  0.05

$means
         X       std   r      LCL      UCL Min Max  Q25 Q50  Q75
1 22.00000  9.137833  21 18.26191 25.73809   4  43 18.0  21 25.0
2 21.37900  8.497583 219 20.22145 22.53654   2  48 16.0  21 27.0
3 17.14815  8.401940  27 13.85146 20.44484   4  32  9.5  15 23.5
4 21.66667 12.971122   9 15.95664 27.37670   2  39 13.0  20 31.0

$comparison
      difference pvalue signif.         LCL      UCL
1 - 2  0.6210046 0.7550          -3.2922091 4.534218
1 - 3  4.8518519 0.0564       .  -0.1322709 9.835975
1 - 4  0.3333333 0.9235          -6.4914578 7.158124
2 - 3  4.2308473 0.0178       *   0.7368444 7.724850
2 - 4 -0.2876712 0.9226          -6.1138493 5.538507
3 - 4 -4.5185185 0.1784         -11.1118932 2.074856

$groups
NULL

attr(,"class")
[1] "group"

Заметим, что p.value позволяет говорить о значимых различиях между группами 2 и 3, с осторожностью – между 1 и 3. Заметим, что с очень высокой уверенностью можно говорить, что различия между соответственно 1 и 2, 1 и 4, 2 и 4 незначимы. Таким образом, люди с неполным высшим образованием (группа “3”) отличаются от группы людей с образованием 8 классов (“1”) или полной средней школы (“2”); группа “4” (люди с высшим образованием) схожа с группами “1” и “2” (неполное и полное среднее образование) и в меньшей степени с группой “3” (неполное высшее) в контексте оценки депресии (bdi).

Для подтверждения этого составим “контрасты” для общей линейной гипотезы и множественных сравнений (glht).

contr <- rbind(
  "1 - 234" = c(-1, 1/3, 1/3, 1/3),
  "2 - 134" = c(1/3, -1, 1/3, 1/3),
  "3 - 124" = c(1/3, 1/3, -1, 1/3),
  "4 - 123" = c(1/3, 1/3, 1/3, -1)
)

GL <- glht(ao, linfct = mcp(group=contr))
summary(GL)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: aov(formula = X ~ group, data = df)

Linear Hypotheses:
             Estimate Std. Error t value Pr(>|t|)  
1 - 234 == 0   -1.935      2.211  -0.875   0.7799  
2 - 134 == 0   -1.107      1.412  -0.785   0.8308  
3 - 124 == 0    4.534      2.044   2.218   0.0932 .
4 - 123 == 0   -1.491      3.027  -0.493   0.9504  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Заметно, как среди прочих сравнений “отличаются” группы “3” и “124”, однако о различии можно говорить лишь с уровнем значимости в 0.1.

Используем разные поправки для множественных сравнений [которые в итоге дадут схожие между собой результаты, подтверждающие выводы выше].

out1 <- LSD.test(ao, "group", p.adj = "bonferroni", group=FALSE)
out2 <- LSD.test(ao, "group", p.adj = "hochberg", group=FALSE)
out3 <- LSD.test(ao, "group", p.adj = "holm", group=FALSE)
out4 <- LSD.test(ao, "group", p.adj = "BH", group=FALSE)
out5 <- LSD.test(ao, "group", p.adj = "fdr", group=FALSE)

out1

$statistics
   MSerror  Df     Mean       CV
  75.70938 272 21.02174 41.39103

$parameters
        test  p.ajusted name.t ntr alpha
  Fisher-LSD bonferroni  group   4  0.05

$means
         X       std   r      LCL      UCL Min Max  Q25 Q50  Q75
1 22.00000  9.137833  21 18.26191 25.73809   4  43 18.0  21 25.0
2 21.37900  8.497583 219 20.22145 22.53654   2  48 16.0  21 27.0
3 17.14815  8.401940  27 13.85146 20.44484   4  32  9.5  15 23.5
4 21.66667 12.971122   9 15.95664 27.37670   2  39 13.0  20 31.0

$comparison
      difference pvalue signif.         LCL       UCL
1 - 2  0.6210046 1.0000          -4.6616672  5.903676
1 - 3  4.8518519 0.3381          -1.8765016 11.580205
1 - 4  0.3333333 1.0000          -8.8798441  9.546511
2 - 3  4.2308473 0.1069          -0.4859078  8.947602
2 - 4 -0.2876712 1.0000          -8.1527635  7.577421
3 - 4 -4.5185185 1.0000         -13.4192935  4.382257

$groups
NULL

attr(,"class")
[1] "group"

pairwise.t.test(data$bdi, data$educat, p.adj = "fdr")

    Pairwise comparisons using t tests with pooled SD 

data:  data$bdi and data$educat 

  1    2    3   
2 0.92 -    -   
3 0.17 0.11 -   
4 0.92 0.92 0.36

P value adjustment method: fdr 

Tukey используется для групп равного объема, поэтому применение здесь некорректно, но ради интереса используем этот метод множественных сравнений средних.

TukeyHSD(ao, "group", ordered=TRUE)

  Tukey multiple comparisons of means
    95% family-wise confidence level
    factor levels have been ordered

Fit: aov(formula = X ~ group, data = df)

$group
         diff        lwr       upr     p adj
2-3 4.2308473 -0.3568649  8.818560 0.0825987
4-3 4.5185185 -4.1387456 13.175783 0.5324422
1-3 4.8518519 -1.6924247 11.396128 0.2235125
4-2 0.2876712 -7.3622447  7.937587 0.9996704
1-2 0.6210046 -4.5171418  5.759151 0.9894140
1-4 0.3333333 -8.6277864  9.294453 0.9996809

---

Непараметрические критерии однородности для независимых выборок

Критерий Манна-Уитни-Вилкоксона (exact=FALSE при объемах выборки больше 30-50)

wilcox.test(X~group, df.2, exact=FALSE, correct=FALSE)

    Wilcoxon rank sum test

data:  X by group
W = 7572, p-value = 0.2038
alternative hypothesis: true location shift is not equal to 0

Критерий Манна-Уитни с поправкой на непрерывность

wilcox.test(X~group, df.2, exact=FALSE, correct=TRUE)

    Wilcoxon rank sum test with continuity correction

data:  X by group
W = 7572, p-value = 0.2041
alternative hypothesis: true location shift is not equal to 0

На основании данного критерия можно сделать вывод, что для данных независимых выборок гипотеза об однородности не отвергается при \(p.value = 0.204\).

Критерий Краскала-Уоллеса для >2 независимых выборок

kruskal.test(X~group, df)

    Kruskal-Wallis rank sum test

data:  X by group
Kruskal-Wallis chi-squared = 5.3869, df = 3, p-value = 0.1456

Данный критерий не отвергает гипотезу об однородности независимых выборок с \(p.value = 0.1456\).

Критерий Краскала с множественными сравнениями и поправкой Бонферрони [результаты позволяют сделать выводы, аналогичные параметрическим методам множественных сравнений].

library(agricolae)
comparison <- with(df, kruskal(X, group, p.adj="bonferroni", group=FALSE, main="HR")); comparison

$statistics
     Chisq Df   p.chisq
  5.386861  3 0.1455644

$parameters
            test  p.ajusted name.t ntr alpha
  Kruskal-Wallis bonferroni  group   4  0.05

$means
         X     rank       std   r Min Max  Q25 Q50  Q75
1 22.00000 144.9286  9.137833  21   4  43 18.0  21 25.0
2 21.37900 141.9658  8.497583 219   2  48 16.0  21 27.0
3 17.14815 104.7593  8.401940  27   4  32  9.5  15 23.5
4 21.66667 140.3889 12.971122   9   2  39 13.0  20 31.0

$comparison
      Difference pvalue Signif.         LCL       UCL
1 - 2   2.962818 1.0000          -45.247856  51.17349
1 - 3  40.169312 0.4994          -21.234930 101.57355
1 - 4   4.539683 1.0000          -79.541538  88.62090
2 - 3  37.206494 0.1342           -5.839515  80.25250
2 - 4   1.576865 1.0000          -70.201473  73.35520
3 - 4 -35.629630 1.0000         -116.859806  45.60055

$groups
NULL

attr(,"class")
[1] "group"

Медианный тест

Median.test(df$X, df$group, correct=TRUE, group=TRUE, console=FALSE)$statistics

---

Анализ зависимых выборок

Данные – dataNF, переменные “BDI.1”, “BDI.2”, “BDI.3” – индекс депрессии в разные моменты времени.

data_dep <- na.omit(dataNF[ , c("BDI.1", "BDI.2", "BDI.3")])

Критерий Стьюдента для зависимых выборок

(paired=TRUE)

t.test(data_dep$BDI.1, data_dep$BDI.2, paired = TRUE)

    Paired t-test

data:  data_dep$BDI.1 and data_dep$BDI.2
t = 11.971, df = 180, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 5.864646 8.179553
sample estimates:
mean of the differences 
               7.022099 

c(t.test(data_dep$BDI.1, data_dep$BDI.2, paired = TRUE)$p.value, t.test(data_dep$BDI.2, data_dep$BDI.3, paired = TRUE)$p.value)

[1] 1.137555e-24 4.269685e-11

По результатам критерия Стьюдента для зависимых выборок гипотеза об однородности изменений во времени отвергается с уровнями значисмости 1.137555e-24 и 4.269685e-11 при сравнении моментов времени 1 и 2; 2 и 3 соответственно (т.е. истинная разность средних не равна 0).

---

Критерий Вилкоксона для зависимых выборок

c(wilcox.test(data_dep$BDI.1, data_dep$BDI.2, paired=TRUE, exact = FALSE)$p.value, wilcox.test(data_dep$BDI.2, data_dep$BDI.3, paired=TRUE, exact = FALSE)$p.value)

[1] 4.195681e-21 1.693720e-11

Аналогично критерию Стьюдента, критерий Вилкоксона отвергает гипотезу об однородности изменений во времени с \(p.value = 4.195681e-21; 1.693720e-11\) при сравненении моментов времени 1 и 2, 2 и 3.

---

Критерий \(\chi^2\)-Фридмана для нескольких выборок

friedman.test(as.matrix(data_dep))

    Friedman rank sum test

data:  as.matrix(data_dep)
Friedman chi-squared = 173.4, df = 2, p-value < 2.2e-16

Непараметрический критерий Фридмана также отвергает гипотезу об однородности “BDI.1”, “BDI.2”, “BDI.3”

---

ANOVA Repeated Measures

Используем модель двухфакторного дисперсионного анализа с повторениями – ANOVA Repeated Measures. В качестве группирующих переменных – “PRCOD.1” и “SEX.1”

dat.AR <- na.omit(dataNF[ , c("PRCOD.1", "SEX.1", "BDI.1", "BDI.3", "BDI.4", "BDI.6")])

k <- 2
m <- ncol(dat.AR)-k; m

[1] 4

dat.AR.T <- data.frame(
    stack(dat.AR[,-seq(k)]),
    sub=as.factor(rep(seq(nrow(dat.AR)), m)),
    gr1=as.factor(rep(dat.AR$PRCOD.1, m)),
    gr2=as.factor(rep(dat.AR$SEX.1, m))
  )

# anova_rm <- aov(values ~ (gr1 + gr2) * ind + Error(sub), dat.AR.T)
anova_rm <- aov(values ~ gr1*gr2 + gr1 + gr2, dat.AR.T)
sarm <- summary(anova_rm); sarm

             Df Sum Sq Mean Sq F value Pr(>F)  
gr1           3    742  247.49   3.545 0.0148 *
gr2           1    108  108.17   1.549 0.2140  
gr1:gr2       3     22    7.33   0.105 0.9571  
Residuals   372  25972   69.82                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Исследуем значимость факторов и эффектов взаимодействия

Df <- sarm[[1]][, 1]
MS <- sarm[[1]][, 3]

p.gr1 <- 1 - pf(MS[1] / MS[4], Df[1], Df[4])
p.gr2 <- 1 - pf(MS[2] / MS[4], Df[2], Df[4])
p.gr12 <- 1 - pf(MS[3] / MS[4], Df[3], Df[4])
c(p.gr1, p.gr2, p.gr12)

[1] 0.01476623 0.21402599 0.95712580

c(sarm[[1]][, 5][1:3])

[1] 0.01476623 0.21402599 0.95712580

\(p > 0.05\), то соответствующий эффект отсутствует.

• Фактор взаимодействия не значим (при нулевой гипотезе об отсутствии эффекта фактора взаимодействия на индекс депрессии BDI \(p.value_{gr1:gr2} = 0.957\), значит, отсутствие эффекта не отвергается, фактор взаимодействия не значим).

• Аналогично отсутствует эффект фактора пола – SEX.1 – на индекс депрессии с \(p.value_{sex.1}= 0.21\), фактор не значим.

• \(p.value_{prcod.1}= 0.0147\), то влияние фактора PRCOD.1 значимо для индекса депрессии, средние в группах по данному фактору значимо различаются (так как нулевая гипотеза гласит об отсутствии влияния).

model.tables(anova_rm, 'mean')

Tables of means
Grand mean
         
11.83421 

 gr1 
    NLTX+Framex NLTX+Placebo Placebo+Framex Placebo+Placebo
          13.89        10.57          11.25            11.4
rep      116.00       136.00          68.00            60.0

 gr2 
    female   male
     12.67  11.52
rep 104.00 276.00

 gr1:gr2 
                 gr2
gr1               female male  
  NLTX+Framex      14.28  13.71
  rep              36.00  80.00
  NLTX+Placebo     12.17  10.22
  rep              24.00 112.00
  Placebo+Framex   12.25  11.04
  rep              12.00  56.00
  Placebo+Placebo  12.09  10.61
  rep              32.00  28.00

Names <- names(table(dat.AR[,1]))
K <- length(Names)

interaction.plot(x.factor=dat.AR.T$ind,
  trace.factor=dat.AR.T$gr1,
  response=dat.AR.T$values,
  fun = mean,
  type = "b", legend = FALSE,
  trace.label ="group",
  xlab = "",
  ylab = 'PRCOD.1',
  lty = seq(K), col = seq(K), pch = 20, lwd = 2
)
legend('topright', Names, lty=seq(K), col=seq(K), cex=0.7, pch=20)

Names <- names(table(dat.AR[,2]))

interaction.plot(x.factor=dat.AR.T$ind,
  trace.factor=dat.AR.T$gr2,
  response=dat.AR.T$values,
  fun = mean,
  type = "b", legend = FALSE,
  trace.label ="group",
  xlab = "",
  ylab = 'SEX.1',
  lty = seq(K), col = seq(K), pch = 20, lwd = 2
)
legend('topright', Names, lty=seq(K), col=seq(K), cex=0.7, pch=20)