Supposons que j'ai mesuré une variable dans les frères et sœurs, qui sont imbriqués dans les familles. La structure des données ressemble à ceci:

valeur de frère et sœur de famille
------ ------- -----
1 1 y_11
1 2 y_12
2 1 y_21
2 2 y_22
2 3 ans_23
... ... ...

Je veux connaître la corrélation entre les mesures prises sur les frères et sœurs au sein de la même famille. La façon habituelle de le faire est de calculer l'ICC sur la base d'un modèle d'interception aléatoire:

res <- lme(yij ~ 1, random = ~ 1 | family, data=dat)
getVarCov(res)[[1]] / (getVarCov(res)[[1]] + res$s^2)

Cela équivaudrait à:

res <- gls(yij ~ 1, correlation = corCompSymm(form = ~ 1 | family), data=dat)

sauf que cette dernière approche permet également un ICC négatif.

Supposons maintenant que j'ai mesuré trois éléments dans des frères et sœurs imbriqués dans des familles. Ainsi, la structure des données ressemble à ceci:

valeur de l'élément de famille de frère
------ ------- ---- -----
1 1 1 y_111
1 1 2 y_112
1 1 3 y_113
1 2 1 y_121
1 2 2 y_122
1 2 3 y_123
2 1 1 y_211
2 1 2 y_212
2 1 3 y_213
2 2 1 a_221
2 2 2 y_222
2 2 3 y_223
2 3 1 a_231
2 3 2 y_232
2 3 3 a_233
... ... ... ...

Maintenant, je veux en savoir plus sur:

1. la corrélation entre les mesures prises sur les frères et sœurs de la même famille pour le même article
2. la corrélation entre les mesures prises sur les frères et sœurs de la même famille pour différents articles

Si je n'avais que des paires de frères et sœurs au sein des familles, je ferais simplement:

res <- gls(yijk ~ item, correlation = corSymm(form = ~ 1 | family), 
           weights = varIdent(form = ~ 1 | item), data=dat)

ce qui me donne une matrice var-cov sur les résidus de la forme:$6 \times 6$

$\left[\begin{array}{ccc|ccc} \sigma^2_1 & \rho_{12} \sigma_1 \sigma_2 & \rho_{13} \sigma_1 \sigma_3 & \phi_{11} \sigma^2_1 & \phi_{12} \sigma_1 \sigma_2 & \phi_{13} \sigma_1 \sigma_3 \\ & \sigma^2_2 & \rho_{23} \sigma_2 \sigma_3 & & \phi_{22} \sigma^2_2 & \phi_{23} \sigma_2 \sigma_3 \\ & & \sigma^2_3 & & & \phi_{33} \sigma^2_3 \\ \hline & & & \sigma^2_1 & \rho_{12} \sigma_1 \sigma_2 & \rho_{13} \sigma_1 \sigma_3 \\ & & & & \sigma^2_2 & \rho_{23} \sigma_2 \sigma_3 \\ & & & & & \sigma^2_3 \\ \end{array}\right]$

$\phi_{jj}$$\phi_{jj'}$

Des idées / suggestions sur la façon dont je pourrais aborder cela? Merci d'avance pour votre aide!

Le package MCMCglmm peut facilement gérer et estimer les structures de covariance et les effets aléatoires. Cependant, il utilise des statistiques bayésiennes qui peuvent intimider les nouveaux utilisateurs. Voir les notes de cours MCMCglmm pour un guide complet sur MCMCglmm, et le chapitre 5 en particulier pour cette question. Je recommande absolument de lire sur l'évaluation de la convergence des modèles et du mélange de chaînes avant d'analyser les données pour de vrai dans MCMCglmm.

library(MCMCglmm)

MCMCglmm utilise des antérieurs, il s'agit d'un précédent wishart inverse non informatif.

p<-list(G=list(
  G1=list(V=diag(2),nu=0.002)),
R=list(V=diag(2),nu=0.002))

Adapter le modèle

m<-MCMCglmm(cbind(x,y)~trait-1,
#trait-1 gives each variable a separate intercept
        random=~us(trait):group,
#the random effect has a separate intercept for each variable but allows and estiamtes the covariance between them.
        rcov=~us(trait):units,
#Allows separate residual variance for each trait and estimates the covariance between them
        family=c("gaussian","gaussian"),prior=p,data=df)

Dans le résumé du modèle, summary(m)la structure G décrit la variance et la covariance des intersections aléatoires. La structure R décrit la variance et la covariance du niveau d'observation de l'interception, qui fonctionnent comme des résidus dans MCMCglmm.

Si vous êtes de persuasion bayésienne, vous pouvez obtenir la distribution postérieure complète des termes de co / variance m$VCV. Notez qu'il s'agit d'écarts après prise en compte des effets fixes.

simuler des données

library(MASS)
n<-3000

#draws from a bivariate distribution
df<-data.frame(mvrnorm(n,mu=c(10,20),#the intercepts of x and y
                   Sigma=matrix(c(10,-3,-3,2),ncol=2)))
#the residual variance covariance of x and y


#assign random effect value
number_of_groups<-100
df$group<-rep(1:number_of_groups,length.out=n)
group_var<-data.frame(mvrnorm(number_of_groups, mu=c(0,0),Sigma=matrix(c(3,2,2,5),ncol=2)))
#the variance covariance matrix of the random effects. c(variance of x,
#covariance of x and y,covariance of x and y, variance of y)

#the variables x and y are the sum of the draws from the bivariate distribution and the random effect
df$x<-df$X1+group_var[df$group,1]
df$y<-df$X2+group_var[df$group,2]

L'estimation de la co / variance d'origine des effets aléatoires nécessite un grand nombre de niveaux pour l'effet aléatoire. Au lieu de cela, votre modèle estimera probablement les co / variances observées qui peuvent être calculées parcov(group_var)

Si vous cherchez à obtenir un "effet familial" et un "effet d'objet", nous pouvons penser qu'il y a des interceptions aléatoires pour ces deux éléments, puis modéliser cela avec le package 'lme4'.

Mais, nous devons d'abord donner à chaque frère et sœur un identifiant unique, plutôt qu'un identifiant unique au sein de la famille.

Ensuite, pour "la corrélation entre les mesures prises sur les frères et sœurs au sein de la même famille pour différents articles", nous pouvons spécifier quelque chose comme:

mod<-lmer(value ~ (1|family)+(1|item), data=family)

Cela nous donnera une interception d'effets fixes pour tous les frères et sœurs, puis deux interceptions d'effets aléatoires (avec variance), pour la famille et l'élément.

Ensuite, pour «la corrélation entre les mesures prises sur les frères et sœurs au sein de la même famille pour le même article», nous pouvons faire la même chose mais simplement sous-définir nos données, nous avons donc quelque chose comme:

mod2<-lmer(value ~ (1|family), data=subset(family,item=="1")) 

Je pense que cela pourrait être une approche plus facile à votre question. Mais, si vous voulez juste l'ICC pour un article ou une famille, le paquet `` psych '' a une fonction ICC () - soyez juste prudent sur la façon dont l'article et la valeur sont fondus dans vos données d'exemple.

Mise à jour

Certains des éléments ci-dessous sont nouveaux pour moi, mais j'ai aimé travailler. Je ne connais vraiment pas l'idée d'une corrélation intraclasse négative. Bien que je vois sur Wikipedia que les «premières définitions ICC» permettaient une corrélation négative avec les données appariées. Mais comme il est le plus couramment utilisé maintenant, l'ICC est compris comme la proportion de la variance totale qui est la variance entre les groupes. Et cette valeur est toujours positive. Bien que Wikipedia ne soit pas la référence la plus fiable, ce résumé correspond à la façon dont j'ai toujours vu ICC utilisé:

Un avantage de ce cadre ANOVA est que différents groupes peuvent avoir différents nombres de valeurs de données, ce qui est difficile à gérer en utilisant les statistiques ICC antérieures. Notez également que cet ICC est toujours non négatif, ce qui permet de l'interpréter comme la proportion de la variance totale qui est «entre les groupes». Cet ICC peut être généralisé pour permettre des effets de covariable, auquel cas l'ICC est interprété comme capturant la similitude intra-classe des valeurs de données ajustées de covariable.

Cela dit, avec des données comme celles que vous avez fournies ici, la corrélation interclasse entre les éléments 1, 2 et 3 pourrait très bien être négative. Et nous pouvons modéliser cela, mais la proportion de la variance expliquée entre les groupes sera toujours positive.

# load our data and lme4
library(lme4)    
## Loading required package: Matrix    

dat<-read.table("http://www.wvbauer.com/fam_sib_item.dat", header=TRUE)

Quel est donc le pourcentage de la variance entre les familles, en contrôlant également la variance entre les groupes entre les groupes d'articles? Nous pouvons utiliser un modèle d'interception aléatoire comme vous l'avez suggéré:

mod<-lmer(yijk ~ (1|family)+(1|item), data=dat)
summary(mod)    
## Linear mixed model fit by REML ['lmerMod']
## Formula: yijk ~ (1 | family) + (1 | item)
##    Data: dat
## 
## REML criterion at convergence: 4392.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6832 -0.6316  0.0015  0.6038  3.9801 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 0.3415   0.5843  
##  item     (Intercept) 0.8767   0.9363  
##  Residual             4.2730   2.0671  
## Number of obs: 1008, groups:  family, 100; item, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)    2.927      0.548   5.342

Nous calculons l'ICC en obtenant la variance des deux intersections à effets aléatoires et des résidus. Nous calculons ensuite le carré de la variance familiale sur la somme des carrés de toutes les variances.

temp<-as.data.frame(VarCorr(mod))$vcov
temp.family<-(temp[1]^2)/(temp[1]^2+temp[2]^2+temp[3]^2)
temp.family    
## [1] 0.006090281

On peut alors faire de même pour les deux autres estimations de variance:

# variance between item-groups
temp.items<-(temp[2]^2)/(temp[1]^2+temp[2]^2+temp[3]^2)
temp.items    
## [1] 0.04015039    
# variance unexplained by groups
temp.resid<-(temp[3]^2)/(temp[1]^2+temp[2]^2+temp[3]^2)
temp.resid    
## [1] 0.9537593    
# clearly then, these will sum to 1
temp.family+temp.items+temp.resid    
## [1] 1

Ces résultats suggèrent que très peu de la variance totale s'explique par la variance entre les familles ou entre les groupes d'articles. Mais, comme indiqué ci-dessus, la corrélation inter-classes entre les éléments pourrait toujours être négative. Commençons par obtenir nos données dans un format plus large:

# not elegant but does the trick
dat2<-cbind(subset(dat,item==1),subset(dat,item==2)[,1],subset(dat,item==3)[,1])
names(dat2)<-c("item1","family","sibling","item","item2","item3")

Nous pouvons maintenant modéliser la corrélation entre, par exemple, item1 et item3 avec une interception aléatoire pour la famille comme auparavant. Mais d'abord, il convient peut-être de rappeler que pour une régression linéaire simple, la racine carrée du r du modèle est la même que le coefficient de corrélation inter-classes (r de Pearson) pour item1 et item2.

# a simple linear regression
mod2<-lm(item1~item3,data=dat2)
# extract pearson's r 
sqrt(summary(mod2)$r.squared)    
## [1] 0.6819125    
# check this 
cor(dat2$item1,dat2$item3)    
## [1] 0.6819125    
# yep, equal

# now, add random intercept to the model
mod3<-lmer(item1 ~ item3 + (1|family), data=dat2)
summary(mod3)    

## Linear mixed model fit by REML ['lmerMod']
## Formula: item1 ~ item3 + (1 | family)
##    Data: dat2
## 
## REML criterion at convergence: 1188.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3148 -0.5348 -0.0136  0.5724  3.2589 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 0.686    0.8283  
##  Residual             1.519    1.2323  
## Number of obs: 336, groups:  family, 100
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -0.07777    0.15277  -0.509
## item3        0.52337    0.02775  18.863
## 
## Correlation of Fixed Effects:
##       (Intr)
## item3 -0.699

La relation entre l'élément1 et l'élément3 est positive. Mais, juste pour vérifier que nous pouvons obtenir une corrélation négative ici, manipulons nos données:

# just going to multiply one column by -1
# to force this cor to be negative

dat2$neg.item3<-dat2$item3*-1
cor(dat2$item1, dat2$neg.item3)    
## [1] -0.6819125    

# now we have a negative relationship
# replace item3 with this manipulated value

mod4<-lmer(item1 ~ neg.item3 + (1|family), data=dat2)
summary(mod4)    

## Linear mixed model fit by REML ['lmerMod']
## Formula: item1 ~ neg.item3 + (1 | family)
##    Data: dat2
## 
## REML criterion at convergence: 1188.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3148 -0.5348 -0.0136  0.5724  3.2589 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 0.686    0.8283  
##  Residual             1.519    1.2323  
## Number of obs: 336, groups:  family, 100
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -0.07777    0.15277  -0.509
## neg.item3   -0.52337    0.02775 -18.863
## 
## Correlation of Fixed Effects:
##           (Intr)
## neg.item3 0.699

Alors oui, la relation entre les articles peut être négative. Mais si nous regardons la proportion de la variance entre les familles dans cette relation, c'est-à-dire ICC (famille), ce nombre sera toujours positif. Comme avant:

temp2<-as.data.frame(VarCorr(mod4))$vcov
(temp2[1]^2)/(temp2[1]^2+temp2[2]^2)    
## [1] 0.1694989

Donc, pour la relation entre l'élément 1 et l'élément 3, environ 17% de cet écart est dû à l'écart entre les familles. Et, nous avons toujours permis qu'il y ait une corrélation négative entre les articles.