Supongamos que he medido alguna variable en hermanos, que están anidados dentro de las familias. La estructura de datos se ve así:

valor de hermano familiar
------ ------- -----
1 1 y_11
1 2 y_12
2 1 y_21
2 2 y_22
2 3 y_23
... ... ...

Quiero saber la correlación entre las medidas tomadas en hermanos dentro de la misma familia. La forma habitual de hacerlo es calcular el ICC basado en un modelo de intercepción aleatoria:

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

Esto sería equivalente a:

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

excepto que el último enfoque también permite una CPI negativa.

Ahora supongamos que he medido tres elementos en hermanos anidados dentro de familias. Entonces, la estructura de datos se ve así:

valor del artículo del hermano de la familia
------ ------- ---- -----
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 y_221
2 2 2 y_222
2 2 3 y_223
2 3 1 y_231
2 3 2 y_232
2 3 3 y_233
... ... ... ...

Ahora, quiero saber sobre:

1. La correlación entre las mediciones tomadas en hermanos dentro de la misma familia para el mismo artículo
2. La correlación entre las mediciones tomadas en hermanos dentro de la misma familia para diferentes artículos

Si solo tuviera pares de hermanos dentro de las familias, simplemente haría:

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

$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'}$

¿Alguna idea / sugerencia de cómo podría abordar esto? ¡Gracias de antemano por cualquier ayuda!

El paquete MCMCglmm puede manejar y estimar fácilmente estructuras de covarianza y efectos aleatorios. Sin embargo, utiliza estadísticas bayesianas que pueden ser intimidantes para los nuevos usuarios. Consulte las Notas del curso MCMCglmm para obtener una guía completa de MCMCglmm, y el capítulo 5 en particular para esta pregunta. Recomiendo leer sobre la evaluación de la convergencia del modelo y la mezcla de cadenas antes de analizar los datos reales en MCMCglmm.

library(MCMCglmm)

MCMCglmm usa priors, este es un previo inverso no informativo.

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

Ajustar el modelo

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)

En el resumen del modelo, summary(m)la estructura G describe la varianza y covarianza de las intersecciones aleatorias. La estructura R describe la varianza del nivel de observación y la covarianza de la intercepción, que funcionan como residuos en MCMCglmm.

Si usted es de persuasión bayesiana, puede obtener la distribución posterior completa de los términos de covarianza m$VCV. Tenga en cuenta que estas son variaciones después de tener en cuenta los efectos fijos.

simular datos

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]

La estimación de la covarianza original de los efectos aleatorios requiere una gran cantidad de niveles para el efecto aleatorio. En cambio, su modelo probablemente estimará las covarianzas observadas que pueden calcularse mediantecov(group_var)

Si está buscando obtener un "efecto familiar" y un "efecto de elemento", podemos pensar que hay interceptaciones aleatorias para ambos y luego modelar esto con el paquete 'lme4'.

Pero, primero tenemos que dar a cada hermano una identificación única, en lugar de una identificación única dentro de la familia.

Luego, para "la correlación entre las mediciones tomadas en hermanos dentro de la misma familia para diferentes elementos", podemos especificar algo como:

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

Esto nos dará una intercepción de efectos fijos para todos los hermanos, y luego dos intercepciones de efectos aleatorios (con variación), para la familia y el elemento.

Luego, para "la correlación entre las mediciones tomadas en hermanos dentro de la misma familia para el mismo elemento", podemos hacer lo mismo pero solo subconjunto nuestros datos, por lo que tenemos algo como:

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

Creo que este podría ser un enfoque más fácil para su pregunta. Pero, si solo desea el ICC para el artículo o la familia, el paquete 'psych' tiene una función ICC (): solo tenga cuidado con la forma en que se fusionan el artículo y el valor en sus datos de ejemplo.

Actualizar

Algunos de los siguientes son nuevos para mí, pero disfruté resolviéndolos. Realmente no estoy familiarizado con la idea de la correlación intraclase negativa. Aunque veo en Wikipedia que "las primeras definiciones de ICC" permitieron una correlación negativa con los datos emparejados. Pero como se usa más comúnmente ahora, ICC se entiende como la proporción de la varianza total que es la varianza entre grupos. Y este valor siempre es positivo. Si bien Wikipedia puede no ser la referencia más autorizada, este resumen corresponde con la forma en que siempre he visto a ICC:

Una ventaja de este marco ANOVA es que diferentes grupos pueden tener diferentes números de valores de datos, lo cual es difícil de manejar utilizando las estadísticas anteriores de ICC. Tenga en cuenta también que este ICC es siempre no negativo, lo que permite interpretarlo como la proporción de la varianza total que es "entre grupos". Este ICC se puede generalizar para permitir efectos covariables, en cuyo caso se interpreta que el ICC captura la similitud dentro de la clase de los valores de datos ajustados por covariable.

Dicho esto, con datos como los que ha proporcionado aquí, la correlación entre clases entre los elementos 1, 2 y 3 podría muy bien ser negativa. Y podemos modelar esto, pero la proporción de la varianza explicada entre los grupos seguirá siendo positiva.

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

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

Entonces, ¿qué porcentaje de la varianza es entre familias, controlando también la varianza entre grupos entre grupos de ítems? Podemos usar un modelo de intercepciones aleatorias como usted sugirió:

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

Calculamos el ICC obteniendo la varianza de las dos intersecciones de efectos aleatorios y de los residuos. Luego calculamos el cuadrado de la varianza familiar sobre la suma de los cuadrados de todas las varianzas.

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

Entonces podemos hacer lo mismo para las otras dos estimaciones de varianza:

# 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

Estos resultados sugieren que muy poca de la varianza total se explica por la varianza entre familias o entre grupos de ítems. Pero, como se señaló anteriormente, la correlación entre clases entre los elementos aún podría ser negativa. Primero, obtengamos nuestros datos en un formato más amplio:

# 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")

Ahora podemos modelar la correlación entre, por ejemplo, el elemento 1 y el elemento 3 con una intercepción aleatoria para la familia como antes. Pero primero, quizás valga la pena recordar que para una regresión lineal simple, la raíz cuadrada del r-cuadrado del modelo es el mismo que el coeficiente de correlación entre clases (r de Pearson) para el elemento 1 y el elemento 2.

# 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 relación es entre item1 y item3 es positiva. Pero, solo para verificar que podamos obtener una correlación negativa aquí, manipulemos nuestros datos:

# 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

Entonces sí, la relación entre los elementos puede ser negativa. Pero si observamos la proporción de variación que hay entre familias en esta relación, es decir, ICC (familia), ese número seguirá siendo positivo. Como antes:

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

Entonces, para la relación entre el elemento 1 y el elemento 3, aproximadamente el 17% de esta variación se debe a la variación entre las familias. Y, todavía hemos permitido que haya una correlación negativa entre los elementos.