Se você deseja obter um "efeito de família" e um "efeito de item", podemos pensar em haver interceptações aleatórias para ambos e modelá-lo com o pacote 'lme4'.
Mas, primeiro, precisamos dar a cada irmão um ID único, em vez de um ID único dentro da família.
Em seguida, para "a correlação entre medidas tomadas em irmãos dentro da mesma família para itens diferentes ", podemos especificar algo como:
mod<-lmer(value ~ (1|family)+(1|item), data=family)
Isso nos dará uma interceptação de efeitos fixos para todos os irmãos e, em seguida, dois efeitos aleatórios interceptam (com variação), para família e item.
Então, para "a correlação entre medições feitas em irmãos dentro da mesma família para o mesmo item", podemos fazer a mesma coisa, mas apenas agrupar nossos dados, para ter algo como:
mod2<-lmer(value ~ (1|family), data=subset(family,item=="1"))
Acho que essa pode ser uma abordagem mais fácil para sua pergunta. Mas, se você deseja apenas o ICC para item ou família, o pacote 'psych' possui uma função ICC () - apenas tenha cuidado com a forma como o item e o valor são derretidos nos dados de exemplo.
Atualizar
Alguns dos itens abaixo são novos para mim, mas gostei de resolver isso. Realmente não estou familiarizado com a ideia de correlação intraclasse negativa. Embora eu veja na Wikipedia que as "primeiras definições da ICC" permitiram uma correlação negativa com os dados emparelhados. Porém, como é mais comumente usado agora, o ICC é entendido como a proporção da variação total que é a variação entre grupos. E esse valor é sempre positivo. Embora a Wikipedia possa não ser a referência mais autorizada, este resumo corresponde a como eu sempre vi o ICC usado:
Uma vantagem dessa estrutura ANOVA é que diferentes grupos podem ter diferentes números de valores de dados, o que é difícil de lidar usando as estatísticas anteriores da ICC. Observe também que este CCI é sempre não negativo, permitindo que seja interpretado como a proporção da variação total "entre grupos". Esse CCI pode ser generalizado para permitir efeitos covariáveis; nesse caso, o CCI é interpretado como capturando o similaridade dentro da classe dos valores de dados ajustados por covariáveis.
Dito isto, com dados como você forneceu aqui, a correlação entre classes entre os itens 1, 2 e 3 pode muito bem ser negativa. E podemos modelar isso, mas a proporção da variação explicada entre os grupos ainda será 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)
Então, qual a porcentagem da variação entre famílias, controlando também a variação entre grupos entre grupos de itens? Podemos usar um modelo de interceptações aleatórias, como você sugeriu:
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 o ICC obtendo a variação dos dois efeitos aleatórios intercepta e dos resíduos. Em seguida, calculamos o quadrado da variação familiar sobre a soma dos quadrados de todas as variações.
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
Podemos então fazer o mesmo para as outras duas estimativas de variação:
# 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
Esses resultados sugerem que muito pouco da variação total é explicado pela variação entre famílias ou entre grupos de itens. Mas, como observado acima, a correlação entre classes entre itens ainda pode ser negativa. Primeiro, vamos obter nossos dados em um formato mais amplo:
# 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")
Agora podemos modelar a correlação entre, por exemplo, item1 e item3 com uma interceptação aleatória para a família como antes. Mas primeiro, talvez valha a pena lembrar que, para uma regressão linear simples, a raiz quadrada do quadrado do r do modelo é igual ao coeficiente de correlação entre classes (r de pearson) para o item1 e o 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
A relação é entre o item1 e o item3 é positivo. Mas, apenas para verificar se podemos obter uma correlação negativa aqui, vamos manipular nossos dados:
# 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
Então, sim, o relacionamento entre os itens pode ser negativo. Mas se observarmos a proporção de variação existente entre as famílias nesse relacionamento, ou seja, ICC (família), esse número ainda será positivo. Como antes:
temp2<-as.data.frame(VarCorr(mod4))$vcov
(temp2[1]^2)/(temp2[1]^2+temp2[2]^2)
## [1] 0.1694989
Portanto, para o relacionamento entre o item1 e o item3, cerca de 17% dessa variação é devido à variação entre famílias. E ainda permitimos que haja uma correlação negativa entre os itens.
var(family)^2/(var(family)^2+var(item)^2)+var(residual)^2)
corresponde? E sim, os ICCs podem ser negativos. Como descrevi no início da minha pergunta, é possível estimar diretamente o ICC com ogls()
modelo, o que permite estimativas negativas. Por outro lado, os modelos de componentes de variação não permitem estimativas negativas.