Załóżmy, że zmierzyłem pewną zmienną w rodzeństwie, które jest zagnieżdżone w rodzinach. Struktura danych wygląda następująco:

wartość rodzeństwa
------ ------- -----
1 1 rok_11
1 2 lata_12
2 1 y_21
2 2 lata_22
2 3 lata_23
... ... ...

Chcę poznać korelację między pomiarami przeprowadzonymi na rodzeństwie w tej samej rodzinie. Typowym sposobem na to jest obliczenie ICC na podstawie modelu losowego przechwytywania:

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

Byłoby to równoważne z:

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

z wyjątkiem tego, że to drugie podejście pozwala również na ujemne ICC.

Załóżmy teraz, że zmierzyłem trzy elementy u rodzeństwa zagnieżdżone w rodzinach. Struktura danych wygląda następująco:

wartość przedmiotu dla rodzeństwa rodzinnego
------ ------- ---- -----
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
... ... ... ...

Teraz chcę się dowiedzieć o:

1. korelacja między pomiarami wykonanymi dla rodzeństwa w tej samej rodzinie dla tego samego elementu
2. korelacja między pomiarami przeprowadzonymi na rodzeństwie w tej samej rodzinie dla różnych przedmiotów

Gdybym miał tylko pary rodzeństwa w rodzinach, po prostu zrobiłbym:

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

co daje mi macierz var-cov na resztach formy:$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]$

na podstawie którego mógłbym łatwo oszacować te korelacje między rodzeństwem (wartości to ICC dla tego samego elementu; wartości to ICC dla różnych elementów). Jednak, jak pokazano powyżej, dla niektórych rodzin mam tylko dwoje rodzeństwa, ale dla innych rodzin więcej niż dwoje. To sprawia, że ​​myślę, że muszę wrócić do modelu typu wariancja-komponenty. Jednak korelacja między elementami może być ujemna, więc nie chcę używać modelu, który ogranicza korelacje jako dodatnie. ϕ j j $\phi_{jj}$$\phi_{jj'}$

Wszelkie pomysły / sugestie, jak mogę do tego podejść? Z góry dziękuję za wszelką pomoc!

Pakiet MCMCglmm może z łatwością obsługiwać i szacować struktury kowariancji i efekty losowe. Wykorzystuje jednak statystyki bayesowskie, które mogą zastraszyć nowych użytkowników. Szczegółowy przewodnik po MCMCglmm znajduje się w notatkach z kursu MCMCglmm, a w szczególności w rozdziale 5 na to pytanie. Absolutnie zalecam przeczytanie o ocenie zbieżności modelu i mieszaniu łańcuchów przed analizowaniem danych rzeczywistych w MCMCglmm.

library(MCMCglmm)

MCMCglmm używa priorów, jest to nieinformacyjny odwrotny życzenie przed.

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

Dopasuj model

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)

W podsumowaniu modelu summary(m)struktura G opisuje wariancję i kowariancję przypadkowych przechwyceń. Struktura R opisuje wariancję poziomu obserwacji i kowariancję punktu przecięcia, które działają jako reszty w MCMCglmm.

Jeśli jesteś perswazją bayesowską, możesz uzyskać całą tylną dystrybucję warunków korelacji / wariancji m$VCV. Zauważ, że są to wariancje po uwzględnieniu ustalonych efektów.

symulować dane

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]

Oszacowanie oryginalnej ko / wariancji efektów losowych wymaga dużej liczby poziomów efektu losowego. Zamiast tego Twój model prawdopodobnie oszacuje zaobserwowane ko / wariancje, które można obliczyćcov(group_var)

Jeśli chcesz uzyskać „efekt rodziny” i „efekt przedmiotu”, możemy pomyśleć o przypadkowych przechwyceniach dla obu z nich, a następnie modelować to za pomocą pakietu „lme4”.

Ale najpierw musimy nadać każdemu rodzeństwu unikalny identyfikator, a nie niepowtarzalny identyfikator w rodzinie.

Następnie dla „korelacji między pomiarami wykonanymi dla rodzeństwa w tej samej rodzinie dla różnych przedmiotów” możemy określić coś takiego:

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

To da nam stałe przechwytywanie efektów dla wszystkich rodzeństwa, a następnie dwa losowe przechwyty efektów (z wariancją) dla rodziny i przedmiotu.

Następnie, w celu „korelacji między pomiarami dokonanymi na rodzeństwie w tej samej rodzinie dla tego samego elementu”, możemy zrobić to samo, ale po prostu podgrupować nasze dane, aby uzyskać coś takiego:

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

Myślę, że może to być łatwiejsze podejście do twojego pytania. Ale jeśli chcesz tylko ICC dla przedmiotu lub rodziny, pakiet „psych” ma funkcję ICC () - zachowaj ostrożność, jeśli chodzi o to, jak przedmiot i wartość są stopione w przykładowych danych.

Aktualizacja

Niektóre z poniższych informacji są dla mnie nowe, ale bardzo mi się podobało. Naprawdę nie znam idei ujemnej korelacji wewnątrzklasowej. Chociaż widzę na Wikipedii, że „wczesne definicje ICC” pozwoliły na ujemną korelację z sparowanymi danymi. Ale jak to jest obecnie najczęściej stosowane, ICC jest rozumiane jako odsetek całkowitej wariancji, która jest wariancją między grupami. Ta wartość jest zawsze dodatnia. Chociaż Wikipedia może nie być najbardziej wiarygodnym odniesieniem, to streszczenie odpowiada temu, w jaki sposób zawsze widziałem ICC:

Zaletą tej struktury ANOVA jest to, że różne grupy mogą mieć różną liczbę wartości danych, co jest trudne do opanowania przy użyciu wcześniejszych statystyk ICC. Należy również zauważyć, że ten ICC jest zawsze nieujemny, co pozwala na interpretację go jako proporcji całkowitej wariancji, która „występuje między grupami”. ICC można uogólnić, aby uwzględnić efekty zmiennych towarzyszących, w którym to przypadku ICC jest interpretowane jako przechwytywanie wewnątrzklasowego podobieństwa wartości danych skorygowanych o zmienne towarzyszące.

To powiedziawszy, przy danych, które tu podałeś, korelacja między klasami między pozycjami 1, 2 i 3 może być bardzo ujemna. Możemy to wymodelować, ale odsetek wariancji wyjaśnionej między grupami nadal będzie dodatni.

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

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

Jaki procent wariancji występuje między rodzinami, kontrolując również wariancję między grupami między grupami przedmiotów? Możemy użyć modelu przechwytywania losowego, tak jak sugerowałeś:

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

Obliczamy ICC, otrzymując wariancję z dwóch przechwyconych efektów losowych i z reszt. Następnie obliczamy kwadrat wariancji rodzinnej na podstawie sumy kwadratów wszystkich wariancji.

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

Możemy wtedy zrobić to samo dla pozostałych dwóch oszacowań wariancji:

# 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

Wyniki te sugerują, że bardzo niewiele całkowitej wariancji można wyjaśnić wariancją między rodzinami lub grupami przedmiotów. Jednak, jak zauważono powyżej, korelacja między klasami między przedmiotami może być nadal ujemna. Najpierw uzyskajmy nasze dane w szerszym formacie:

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

Teraz możemy modelować korelację między, na przykład, item1 i item3 za pomocą losowego przechwytywania dla rodziny, jak poprzednio. Ale po pierwsze, być może warto przypomnieć, że dla prostej regresji liniowej pierwiastek kwadratowy r-kwadrat modelu jest taki sam, jak współczynnik korelacji między klasami (r Pearsona dla item1 i 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

Zależność między pozycją 1 a pozycją 3 jest dodatnia. Aby jednak sprawdzić, czy możemy uzyskać ujemną korelację, zmanipulujmy nasze dane:

# 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

Tak więc związek między przedmiotami może być ujemny. Ale jeśli spojrzymy na proporcję wariancji między rodzinami w tym związku, tj. ICC (rodzina), liczba ta nadal będzie dodatnia. Jak wcześniej:

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

Tak więc dla relacji między item1 i item3 około 17% tej wariancji wynika z wariancji między rodzinami. Nadal pozwalaliśmy na istnienie ujemnej korelacji między przedmiotami.