Как сделать: корреляция с «блоками» (или - «повторными мерами»?!)?

Question

У меня есть следующие настройки для анализа: У нас есть около 150 предметов, и для каждого предмета мы провели пару тестов (в разных условиях) 18 раз. 18 различных условий теста являются взаимодополняющими, таким образом, что если мы будем усреднять результаты тестов (для каждого субъекта), мы не получим никакой корреляции между тестами (между субъектами). То, что мы хотим знать, это корреляция (и значение P) между тестами в рамках предметов, но по всем предметам.

Способ, которым я сделал это сейчас, состоял в том, чтобы выполнить корреляцию для каждого субъекта, а затем посмотреть на распределение полученных корреляций, чтобы увидеть, отличается ли это среднее от 0. Но я подозреваю, что может быть лучший способ ответить на тот же вопрос (кто-то сказал мне что-то о «географической корреляции», но поверхностный поиск не помог).

p.s: Я понимаю, что здесь может быть место для создания какой-то смешанной модели, но я бы предпочел представить «корреляцию», и я не уверен, как извлечь такой вывод из смешанной модели.

Кроме того, вот короткий фиктивный код, чтобы дать представление о том, о чем я говорю:

attach(longley)
N <- length(Unemployed)
block <- c(
        rep( "a", N),
        rep( "b", N),
        rep( "c", N)
        )

Unemployed.3 <- c(Unemployed + rnorm(1),
                    Unemployed + rnorm(1),
                    Unemployed + rnorm(1))

GNP.deflator.3 <- c(GNP.deflator + rnorm(1),
                    GNP.deflator + rnorm(1),
                    GNP.deflator + rnorm(1))

cor(Unemployed, GNP.deflator)
cor(Unemployed.3, GNP.deflator.3)
cor(Unemployed.3[block == "a"], GNP.deflator.3[block == "a"])
cor(Unemployed.3[block == "b"], GNP.deflator.3[block == "b"])
cor(Unemployed.3[block == "c"], GNP.deflator.3[block == "c"])
(I would like to somehow combine the last three correlations...)

Любые идеи будут приветствоваться.

Лучший, Tal

Answer 1

Я согласен с Тристаном - вы ищете ICC. Единственное отличие от стандартных реализаций состоит в том, что два оценщика (тесты) оценивают каждого субъекта повторно. Там может быть реализация, которая позволяет это. Между тем вот еще один подход для получения корреляции.

Вы можете использовать «общие линейные модели», которые являются обобщениями линейных моделей, которые явно допускают корреляцию между невязками. Код ниже реализует это с помощью функции gls пакета nlme. Я уверен, что есть и другие способы. Чтобы использовать эту функцию, мы должны сначала преобразовать данные в «длинный» формат. Я также изменил имена переменных на x и y для простоты. Я также использовал +rnorm(N) вместо +rnorm(1) в вашем коде, потому что, я думаю, вы имели в виду.

library(reshape)
library(nlme)
dd <- data.frame(x=Unemployed.3, y=GNP.deflator.3, block=factor(block))
dd$occasion <- factor(rep(1:N, 3))  # variable denoting measurement occasions
dd2 <- melt(dd, id=c("block","occasion"))  # reshape

# fit model with the values within a measurement occasion correlated
#   and different variances allowed for the two variables
mod <- gls(value ~ variable + block, data=dd2, 
           cor=corSymm(form=~1|block/occasion), 
           weights=varIdent(form=~1|variable))  
# extract correlation
mod$modelStruct$corStruct

В платформе моделирования вы можете использовать критерий отношения правдоподобия, чтобы получить значение p. nlme также может дать вам доверительный интервал:

mod2 <- gls(value ~ variable + block, data=dd2, 
           weights=varIdent(form=~1|variable))  
anova(mod, mod2)   # likelihood-ratio test for corr=0

intervals(mod)$corStruct  # confidence interval for the correlation

Answer 2

Если я правильно понимаю ваш вопрос, вы заинтересованы в вычислении внутриклассовой корреляции между несколькими тестами. В пакете psy есть реализация, хотя я ее не использовал.

Если вы хотите сделать вывод относительно оценки корреляции, вы можете начать тестирование с предметов. Просто убедитесь, что тесты хранятся вместе для каждого образца.

Answer 3

Я не эксперт, но это выглядит для меня как то, что вы хотите. Он автоматизирован, имеет короткий код, дает те же корреляции, что и в приведенном выше примере, и выдает p-значения.

> df = data.frame(block=block, Unemployed=Unemployed.3,
+ GNP.deflator=GNP.deflator.3)
> require(plyr)
Loading required package: plyr
> ddply(df, "block", function(x){
+   as.data.frame(
+     with(x,cor.test(Unemployed, GNP.deflator))[c("p.value","estimate")]
+ )})
  block    p.value  estimate
1     a 0.01030636 0.6206334
2     b 0.01030636 0.6206334
3     c 0.01030636 0.6206334

Чтобы увидеть все детали, сделайте следующее:

> dlply(df, "block", function(x){with(x,cor.test(Unemployed, GNP.deflator))})
$a

    Pearson's product-moment correlation

data:  Unemployed and GNP.deflator 
t = 2.9616, df = 14, p-value = 0.01031
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 0.1804410 0.8536976 
sample estimates:
      cor 
0.6206334 


$b

    Pearson's product-moment correlation

data:  Unemployed and GNP.deflator 
t = 2.9616, df = 14, p-value = 0.01031
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 0.1804410 0.8536976 
sample estimates:
      cor 
0.6206334 


$c

    Pearson's product-moment correlation

data:  Unemployed and GNP.deflator 
t = 2.9616, df = 14, p-value = 0.01031
alternative hypothesis: true correlation is not equal to 0 
95 percent confidence interval:
 0.1804410 0.8536976 
sample estimates:
      cor 
0.6206334 


attr(,"split_type")
[1] "data.frame"
attr(,"split_labels")
  block
1     a
2     b
3     c