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

Question

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

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

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

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

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...)

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

Лучший, Тал

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

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

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

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