Come fare: Correlazione con “blocchi” (o - “misure ripetute”?!)?

Question

Ho la seguente messa a punto da analizzare: Abbiamo circa 150 soggetti, e per ogni soggetto è stata eseguita una coppia di test (in condizioni diverse) 18 volte. Le 18 diverse condizioni della prova sono complementari, in modo tale che se dove media nel test (per ciascun soggetto), si otterrebbe alcuna correlazione tra le prove (tra soggetti). Ciò che vogliamo sapere è la correlazione (e il valore P) tra le prove, in entro i soggetti, ma su tutti gli argomenti.

Il modo in cui ho fatto questo, ormai è stato quello di eseguire la correlazione per ogni soggetto, e poi guardare la distribuzione delle correlazioni ricevuto così per vedere se si tratta di media è diverso quindi 0. Ma ho il sospetto che ci potrebbe essere un modo migliore per rispondere alla stessa domanda (qualcuno mi ha detto qualcosa su "correlazione geografica", ma una ricerca superficiale non ha aiutato).

ps:. Capisco che ci potrebbe essere un posto qui per fare una sorta di modello misto, ma io preferirei di presentare una "correlazione", e non sono sicuro di come estrarre una tale uscita da un modello misto

Inoltre, ecco un codice fittizio breve per dare un'idea di cosa sto parlando:

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

Tutte le idee saranno accolti.

Best, Tal

Answer 1

Sono d'accordo con Tristan - siete alla ricerca di ICC. L'unica differenza rispetto implementazioni standard è che i due valutatori (test) valutare ciascun soggetto ripetutamente. Ci potrebbe essere un'implementazione che permette questo. Nel frattempo ecco un altro approccio per ottenere la correlazione.

È possibile utilizzare "modelli lineari generali", che sono generalizzazioni di modelli lineari che consentono esplicitamente la correlazione tra i residui. Il seguente codice implementa questa utilizzando la funzione gls del pacchetto nlme. Sono sicuro che ci sono altri modi. Per utilizzare questa funzione dobbiamo rimodellare prima i dati in un formato "lungo". Ho anche cambiato i nomi delle variabili per x e y per semplicità. Ho anche usato +rnorm(N) invece di +rnorm(1) nel codice, perché è quello che penso che volevi dire.

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

Nel quadro di modellazione è possibile utilizzare un test di rapporto di verosimiglianza per ottenere un p-value. nlme può anche dare un intervallo di confidenza:

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

Se ho capito bene la tua domanda, siete interessati a calcolare le href="http://en.wikipedia.org/wiki/Intraclass_correlation" correlazione intraclasse tra multipla test. C'è un'implementazione nel pacchetto psy , anche se ho non utilizzato.

Se si desidera eseguire l'inferenza sulla stima di correlazione, è possibile il bootstrap dei soggetti. Basta fare in modo di tenere insieme i test per ogni campione.

Answer 3

Non sono un esperto, ma questo mi sembra quello che vuoi. E 'automatizzato, breve per il codice, dà gli stessi correlazioni come vostro esempio di cui sopra, e produce valori di 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

Per visualizzare tutti i dettagli, fare questo:

> 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