Cómo hacer: Correlación con “bloques” (o “-” medidas repetidas?!)?

Question

I tiene la siguiente configuración para analizar: Tenemos alrededor de 150 sujetos, y para cada sujeto se realizó un par de pruebas (en condiciones diferentes) 18 veces. Las 18 condiciones diferentes de la prueba son complementarias, de tal manera de modo que si dónde promedio durante las pruebas (para cada sujeto), que se pueden conseguir ninguna correlación entre las pruebas (entre los sujetos). Lo que queremos saber es la correlación (y el valor P) entre las pruebas, en el plazo de los sujetos, pero sobre todo los temas.

La forma en que hice esto por ahora ha sido realizar la correlación para cada sujeto, y luego ver la distribución de las correlaciones recibido hasta para ver si es media es diferente a 0. Pero sospecho que puede haber una mejor manera de responder a la misma pregunta (alguien me dijo algo acerca de "correlación geográfica", pero una búsqueda superficial no ayudó).

ps:. Entiendo que puede haber un lugar para hacer algún tipo de modelo mixto, pero yo preferiría presentar una "correlación", y no estoy seguro de cómo extraer dicha salida de un modelo mixto

Además, aquí es un código ficticio corta para dar una idea de lo que estoy hablando:

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

será recibido ninguna idea.

Best, Tal

Answer 1

Estoy de acuerdo con Tristán - que busca la CPI. La única diferencia con implementaciones estándar es que los dos evaluadores (pruebas) evalúan cada sujeto repetidamente. Es posible que haya una aplicación que permite eso. Mientras tanto aquí es otro enfoque para obtener la correlación.

Se puede usar "modelos lineales generales", que son generalizaciones de modelos lineales que permiten explícitamente correlación entre los residuales. El código siguiente implementa esta función utilizando el gls del paquete nlme. Estoy seguro de que hay otras maneras también. Para utilizar esta función hay que formar de nuevo la primera de datos en un formato de "largo". También ha cambiado los nombres de las variables a x y y por simplicidad. También utilicé +rnorm(N) en lugar de +rnorm(1) en su código, porque eso es lo que creo que quería decir.

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

En el marco de modelado se puede utilizar una prueba de razón de verosimilitud para obtener un valor de p. nlme también puede darle un intervalo de confianza:

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

Si entiendo bien su pregunta, usted está interesado en el cálculo de la correlación intraclase entre múltiples pruebas. Hay una aplicación en el paquete psy , aunque no tengo no utilizado.

Si desea realizar inferencia en la estimación de correlación, que podría arrancar los sujetos. Sólo asegúrese de mantener juntos las pruebas para cada muestra.

Answer 3

No soy un experto, pero esto me parece lo que quiero. Se automatizado, corta para código, da las mismas correlaciones como su ejemplo de arriba, y produce los valores de 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

Para ver todos los detalles, hacer esto:

> 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