Предположим, я измерил некоторую переменную у братьев и сестер, которые вложены в семьи. Структура данных выглядит следующим образом:

семейный брат
------ ------- -----
1 1 год_11
1 2 год_12
2 1 год_21
2 2 y_22
2 3 y_23
... ... ...

Я хочу знать корреляцию между измерениями, сделанными на братьях и сестрах в пределах одной семьи. Обычный способ сделать это - рассчитать ICC на основе модели случайного перехвата:

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

Это будет эквивалентно:

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

за исключением того, что последний подход также допускает отрицательный ICC.

Теперь предположим, что я измерил три элемента в братьях и сестрах, вложенных в семьи. Итак, структура данных выглядит так:

стоимость предмета родного брата
------ ------- ---- -----
1 1 1 год_111
1 1 2 y_112
1 1 3 y_113
1 2 1 год_121
1 2 2 y_122
1 2 3 y_123
2 1 1 год_211
2 1 2 y_212
2 1 3 y_213
2 2 1 год_221
2 2 2 y_222
2 2 3 y_223
2 3 1 год_231
2 3 2 y_232
2 3 3 y_233
... ... ... ...

Теперь я хочу узнать о:

1. корреляция между измерениями, взятыми на братьев и сестер в пределах одной семьи для одного и того же элемента
2. корреляция между измерениями, взятыми на братьев и сестер в пределах одного и того же семейства для разных предметов

Если бы в семьях были только пары братьев и сестер, я бы просто сделал:

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

Любые идеи / предложения о том, как я мог бы подойти к этому? Заранее благодарю за любую помощь!

Пакет MCMCglmm может легко обрабатывать и оценивать ковариационные структуры и случайные эффекты. Однако он использует байесовскую статистику, которая может быть пугающей для новых пользователей. См. Примечания к курсу MCMCglmm для подробного руководства по MCMCglmm, и главу 5, в частности, для этого вопроса. Я настоятельно рекомендую ознакомиться с оценкой сходимости моделей и смешивания цепей, прежде чем анализировать данные в реальном времени в MCMCglmm.

library(MCMCglmm)

MCMCglmm использует априоры, это неинформативное обратное пожелание ранее.

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

Подходит модель

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)

В сводке модели summary(m)структура G описывает дисперсию и ковариацию случайных перехватов. Структура R описывает дисперсию уровня наблюдения и ковариацию перехвата, которые функционируют как невязки в MCMCglmm.

Если у вас есть байесовское убеждение, вы можете получить полное заднее распределение слагаемых / дисперсии m$VCV. Обратите внимание, что это отклонения после учета фиксированных эффектов.

симулировать данные

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]

Для оценки исходного коэффициента случайных эффектов требуется большое количество уровней случайного эффекта. Вместо этого ваша модель, вероятно, будет оценивать наблюдаемые ко / дисперсии, которые можно рассчитатьcov(group_var)

Если вы хотите получить «семейный эффект» и «эффект предмета», мы можем подумать о случайных перехватах для них обоих, а затем смоделировать это с помощью пакета «lme4».

Но сначала мы должны дать каждому брату уникальный идентификатор, а не уникальный идентификатор в семье.

Затем для «корреляции между измерениями, проведенными на братьях и сестрах в пределах одного и того же семейства для разных предметов», мы можем указать что-то вроде:

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

Это даст нам фиксированный перехват эффектов для всех братьев и сестер, а затем два случайных перехвата эффектов (с дисперсией) для семьи и предмета.

Затем, для «корреляции между измерениями, проведенными для братьев и сестер в пределах одного и того же семейства для одного и того же предмета», мы можем сделать то же самое, но просто подналожить наши данные, поэтому у нас есть что-то вроде

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

Я думаю, что это может быть более простым подходом к вашему вопросу. Но если вы просто хотите получить ICC для предмета или семейства, в пакете 'psych' есть функция ICC () - просто будьте осторожны с тем, как элемент и ценность расплавляются в данных вашего примера.

Обновить

Некоторое из нижеприведенного является новым для меня, но мне понравилось работать над ним. Я действительно не знаком с идеей отрицательной внутриклассовой корреляции. Хотя я вижу в Википедии, что «ранние определения ICC» допускали отрицательную корреляцию с парными данными. Но так как он чаще всего используется сейчас, ICC понимается как доля от общей дисперсии, которая является дисперсией между группами. И это значение всегда положительно. Хотя Википедия, возможно, не самая авторитетная ссылка, это резюме соответствует тому, как я всегда видел использование ICC:

Преимущество этой структуры ANOVA состоит в том, что разные группы могут иметь разное количество значений данных, что трудно обрабатывать с использованием более ранней статистики ICC. Также обратите внимание, что этот ICC всегда неотрицателен, что позволяет интерпретировать его как долю общей дисперсии, которая «между группами». Этот ICC может быть обобщен для учета ковариатных эффектов, и в этом случае ICC интерпретируется как захват сходства внутри класса значений данных, скорректированных по ковариату.

Тем не менее, с данными, которые вы привели здесь, межклассовая корреляция между пунктами 1, 2 и 3 вполне может быть отрицательной. И мы можем смоделировать это, но доля дисперсии, объясненной между группами, все еще будет положительной.

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

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

Так какой процент различий между семьями, учитывая также дисперсию между группами между группами предметов? Мы можем использовать модель случайных перехватов, как вы предложили:

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

Мы рассчитываем ICC, получая дисперсию от двух случайных перехватов эффектов и от остатков. Затем мы вычисляем квадрат семейной дисперсии по сумме квадратов всех дисперсий.

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

Затем мы можем сделать то же самое для двух других оценок дисперсии:

# 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

Эти результаты показывают, что очень мало общей дисперсии объясняется разницей между семействами или группами предметов. Но, как отмечалось выше, межклассовая корреляция между предметами все еще может быть отрицательной. Сначала давайте получим наши данные в более широком формате:

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

Теперь мы можем смоделировать корреляцию между, например, item1 и item3 со случайным перехватом для family, как и раньше. Но сначала, возможно, стоит напомнить, что для простой линейной регрессии квадратный корень из r-квадрата модели совпадает с коэффициентом корреляции между классами (r Пирсона) для item1 и 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

Связь между item1 и item3 является положительной. Но, просто чтобы проверить, что мы можем получить здесь отрицательную корреляцию, давайте манипулируем нашими данными:

# 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

Так что да, отношения между предметами могут быть отрицательными. Но если мы посмотрим на долю дисперсии между семьями в этих отношениях, то есть ICC (семья), это число все равно будет положительным. Как прежде:

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

Таким образом, для отношений между item1 и item3, около 17% этой разницы происходит из-за различий между семействами. И мы по-прежнему допускаем наличие отрицательной корреляции между предметами.