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.
var(family)^2/(var(family)^2+var(item)^2)+var(residual)^2)
odpowiada? I tak, ICC mogą być ujemne. Jak opisałem na początku mojego pytania, można bezpośrednio oszacować ICC za pomocągls()
modelu, co pozwala na oszacowania ujemne. Z drugiej strony modele składowe wariancji nie pozwalają na ujemne oszacowania.