Comment faire: corrélation avec les « blocs » (ou « - mesures répétées »?!)?

StackOverflow https://stackoverflow.com/questions/2336056

  •  22-09-2019
  •  | 
  •  

Question

J'ai la configuration suivante pour analyser: Nous avons environ 150 sujets, et pour chaque sujet nous avons effectué deux essais (dans des conditions différentes) 18 fois. Les 18 différentes conditions de l'essai sont complémentaires, de telle manière de sorte que si nous où en moyenne au cours des essais (pour chaque sujet), nous obtiendrions pas de corrélation entre les essais (entre les sujets). Ce que nous voulons savoir est la corrélation (et la valeur P) entre les essais, à l'intérieur des sujets, mais sur tous les sujets.

La façon dont je l'ai fait maintenant était d'effectuer la corrélation pour chaque sujet, puis examiner la répartition des corrélations reçues afin de voir si elle est moyenne est différente de 0. Mais je pense qu'il pourrait y avoir un meilleur moyen pour répondre à la même question (quelqu'un m'a dit quelque chose sur « corrélation géographique », mais une recherche superficielle n'a pas aidé).

ps. Je comprends qu'il pourrait y avoir une place ici pour faire une sorte de modèle mixte, mais je préfère présenter une « corrélation », et ne suis pas sûr de savoir comment extraire une telle sortie d'un modèle mixte

En outre, voici un petit code factice pour donner une idée de ce dont je parle: 

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

Toutes les idées seront les bienvenus.

Best, Tal

Était-ce utile?

La solution

Je suis d'accord avec Tristan - vous cherchez la CPI. La seule différence par rapport à des implémentations standard est que les deux évaluateurs (tests) évaluent chaque sujet à plusieurs reprises. Il pourrait y avoir une implémentation qui permet. En attendant, est ici une autre approche pour obtenir la corrélation.

Vous pouvez utiliser des « modèles linéaires généraux », qui sont des généralisations des modèles linéaires qui permettent explicitement la corrélation entre les résidus. Le code ci-dessous met en œuvre grâce à la fonction de gls du paquet nlme. Je suis sûr qu'il ya d'autres façons. Pour utiliser cette fonction, nous devons d'abord remodeler les données dans un format « long ». J'ai aussi changé les noms de variables à x et y pour la simplicité. J'ai aussi utilisé +rnorm(N) au lieu de +rnorm(1) dans votre code, parce que ce que je pense que vous vouliez 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

Dans le cadre de la modélisation, vous pouvez utiliser un test de rapport de vraisemblance pour obtenir une valeur de p. nlme peut aussi vous donner un intervalle de confiance: 

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

Autres conseils

Je ne suis pas expert, mais cela me semble que ce que vous voulez. Il est automatisé, court au code, donne les mêmes corrélations que votre exemple ci-dessus, et produit des valeurs 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

Pour voir tous les détails, faites ceci: 

> 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
Non affilié à StackOverflow
scroll top