如何操作：與「塊」（或 - 「重複測量」？！）相關？

Question

我有以下設置來分析： 我們有大約150個科目，並且對於每個科目我們都進行了18次測試（在不同條件下）。 18個不同的測試條件是互補的，這樣，如果我們在測試中的平均值（對於每個對象），測試之間（對象之間）沒有相關性。 我們希望知道的是測試之間的相關性（和P值），在受試者內部，但是在所有受試者中。

我現在這樣做的方式是執行爲每個主題的相關性，再看看相關的分佈收到這麼看它是否意味着是不同的，那麼0 但我懷疑有可能是一個更好的回答同樣問題的方式（有人對我說了一些關於「地理相關性」的內容，但淺搜索沒有幫助）。 p：我知道這裏可能有一個地方可以做某種混合模型，但我更願意提出一個「相關性」，我不知道如何從混合模型中提取這樣的輸出。

而且，這裏是一個短路僞代碼，給什麼我談論的一個想法：

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。與標準實現唯一的區別是兩個評估者（測試）反覆評估每個主題。可能有一個允許的實現。與此同時，這裏是另一種獲得相關性的方法。

您可以使用「通用線性模型」，它是線性模型的推廣，明確允許殘差之間的相關性。下面的代碼使用nlme包的gls函數來實現此功能。我相信還有其他的方法。要使用此功能，我們必須首先將數據重塑爲「長」格式。爲了簡單起見，我還將變量名稱更改爲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

我不是專家，但在我看來像你想要什麼。它是自動的，代碼簡潔，與上面的示例給出了相同的相關性，並生成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 

Answer 3

如果我理解你的問題正確，你有興趣在計算多個測試之間的intraclass correlation。儘管我沒有使用它，但在psy包中有一個實現。

如果你想對關聯估計進行推理，你可以引導主體。只要確保爲每個樣品保留測試。