ICC jako oczekiwana korelacja między dwiema losowo losowanymi jednostkami, które są w tej samej grupie

12

W modelowaniu wielopoziomowym korelacja wewnątrzklasowa jest często obliczana na podstawie analizy ANOVA z efektami losowymi

yjajot=γ00+ujot+mijajot

gdzie to poziomu 2, a to reszty poziomu 1. Następnie otrzymujemy oszacowania, i dla wariancji odpowiednio i , i je do następującego równania:ujotmijajotσ^u2)σ^mi2)ujotmijajot

ρ=σ^u2)σ^u2)+σ^mi2)

Hox (2002) pisze na s. 15, że

Korelację wewnątrzklasową ρ można również interpretować jako oczekiwaną korelację między dwiema losowymi jednostkami, które są w tej samej grupie

Jest tutaj pytanie , które zadaje zaawansowane pytanie (dlaczego jest dokładnie takie samo, a nie w przybliżeniu równe) i otrzymuje zaawansowaną odpowiedź.

Chciałbym jednak zadać o wiele prostsze pytanie.

Pytanie: Co to w ogóle oznacza mówienie o korelacji między dwiema losowymi jednostkami, które są w tej samej grupie?

Rozumiem, że korelacja wewnątrzklasowa działa na grupy, a nie na sparowane dane. Jednak nadal nie rozumiem, jak można obliczyć korelację, gdybyśmy mieli tylko dwie losowo losowane jednostki z tej samej grupy. Na przykład, jeśli spojrzę na wykresy punktowe na stronie Wikipedii dla ICC , mamy wiele grup i wiele punktów w każdej grupie.

user1205901 - Przywróć Monikę
źródło

Odpowiedzi:

10

Najłatwiej jest zobaczyć równoważność, jeśli weźmie się pod uwagę przypadek, w którym są tylko dwie osoby na grupę. Przejdźmy więc do konkretnego przykładu (użyję do tego R):

dat <- read.table(header=TRUE, text = "
group person   y
1     1        5
1     2        6
2     1        3
2     2        2
3     1        7
3     2        9
4     1        2
4     2        2
5     1        3
5     2        5
6     1        6
6     2        9
7     1        4
7     2        2
8     1        8
8     2        7")

Mamy więc 8 grup po 2 osoby. Teraz dopasujmy model ANOVA z efektami losowymi:

library(nlme)
res <- lme(y ~ 1, random = ~ 1 | group, data=dat, method="ML")

Na koniec obliczmy ICC:

getVarCov(res)[1] / (getVarCov(res)[1] + res$sigma^2)

Daje to: 0.7500003(dokładnie 0,75, ale w procedurze szacowania istnieje pewne niewielkie wrażenie numeryczne).

Teraz przekształćmy dane z długiego formatu na szeroki format:

dat <- as.matrix(reshape(dat, direction="wide", v.names="y", idvar="group", timevar="person"))

Teraz wygląda to tak:

   group y.1 y.2
1      1   5   6
3      2   3   2
5      3   7   9
7      4   2   2
9      5   3   5
11     6   6   9
13     7   4   2
15     8   8   7

A teraz obliczyć korelację między y.1i y.2:

cor(dat[,2], dat[,3])

Daje to: 0.8161138

Czekaj, co? Co tu się dzieje? Czy nie powinno to być 0,75? Nie do końca! To, co obliczyłem powyżej, to nie ICC ( wewnątrzklasowy współczynnik korelacji), ale zwykły współczynnik korelacji iloczynu Pearsona, który jest współczynnikiem korelacji międzyklasowej . Zauważ, że w danych o długim formacie jest całkowicie arbitralne, kto jest osobą 1, a kto jest osobą 2 - pary są nieuporządkowane. Możesz przetasować dane w grupach i uzyskasz te same wyniki. Ale w danych wielkoformatowych nie jest arbitralne, kto jest wymieniony poniżej, y.1a kto wymieniony poniżej y.2. Gdybyś przełączył niektóre osoby, uzyskałbyś inną korelację (z wyjątkiem gdybyś przerzucił je wszystkie - to odpowiada tocor(dat[,3], dat[,2])co oczywiście nadal daje 0.8161138).

Fisher wskazał na małą sztuczkę, aby uzyskać ICC z danymi szerokoformatowymi. Niech każda para zostanie uwzględniona dwukrotnie, w obu zamówieniach, a następnie obliczyć korelację:

dat <- rbind(dat, dat[,c(1,3,2)])
cor(dat[,2], dat[,3])

Daje to: 0.75.

Jak więc widać, ICC jest tak naprawdę współczynnikiem korelacji - dla „niesparowanych” danych dwóch osób z tej samej grupy.

Jeśli w grupie są więcej niż dwie osoby, nadal możesz myśleć o ICC w ten sposób, z tym wyjątkiem, że istnieje więcej sposobów tworzenia par osób w grupach. ICC jest wówczas korelacją między wszystkimi możliwymi parami (ponownie w nieuporządkowany sposób).

Wolfgang
źródło
7

y

Najpierw ładujemy zestaw danych @ Wolfganga (nie pokazano tutaj). Teraz zdefiniujmy prostą funkcję R, która pobiera ramkę danych i zwraca pojedynczą losowo wybraną parę obserwacji z tej samej grupy:

get_random_pair <- function(df){
  # select a random row
  i <- sample(nrow(df), 1)
  # select a random other row from the same group
  # (the call to rep() here is admittedly odd, but it's to avoid unwanted
  # behavior when the first argument to sample() has length 1)
  j <- sample(rep(setdiff(which(dat$group==dat[i,"group"]), i), 2), 1)
  # return the pair of y-values
  c(df[i,"y"], df[j,"y"])
}

Oto przykład tego, co otrzymamy, jeśli wywołamy tę funkcję 10 razy w zestawie danych @ Wolfgang:

test <- replicate(10, get_random_pair(dat))
t(test)
#       [,1] [,2]
#  [1,]    9    6
#  [2,]    2    2
#  [3,]    2    4
#  [4,]    3    5
#  [5,]    3    2
#  [6,]    2    4
#  [7,]    7    9
#  [8,]    5    3
#  [9,]    5    3
# [10,]    3    2

Aby oszacować ICC, wywołujemy tę funkcję wiele razy, a następnie obliczamy korelację między dwiema kolumnami.

random_pairs <- replicate(100000, get_random_pair(dat))
cor(t(random_pairs))
#           [,1]      [,2]
# [1,] 1.0000000 0.7493072
# [2,] 0.7493072 1.0000000

Tę samą procedurę można zastosować, bez żadnych modyfikacji, do zbiorów danych z grupami dowolnej wielkości. Na przykład stwórzmy zestaw danych składający się ze 100 grup po 100 obserwacji każda, z prawdziwym ICC ustawionym na 0,75, jak w przykładzie @ Wolfganga.

set.seed(12345)
group_effects <- scale(rnorm(100))*sqrt(4.5)
errors <- scale(rnorm(100*100))*sqrt(1.5)
dat <- data.frame(group = rep(1:100, each=100),
                  person = rep(1:100, times=100),
                  y = rep(group_effects, each=100) + errors)

stripchart(y ~ group, data=dat, pch=20, col=rgb(0,0,0,.1), ylab="group")

wprowadź opis zdjęcia tutaj

Szacując ICC na podstawie składników wariancji z modelu mieszanego, otrzymujemy:

library("lme4")
mod <- lmer(y ~ 1 + (1|group), data=dat, REML=FALSE)
summary(mod)
# Random effects:
#  Groups   Name        Variance Std.Dev.
#  group    (Intercept) 4.502    2.122   
#  Residual             1.497    1.223   
# Number of obs: 10000, groups:  group, 100

4.502/(4.502 + 1.497)
# 0.7504584

A jeśli zastosujemy losową procedurę parowania, otrzymamy

random_pairs <- replicate(100000, get_random_pair(dat))
cor(t(random_pairs))
#           [,1]      [,2]
# [1,] 1.0000000 0.7503004
# [2,] 0.7503004 1.0000000

co ściśle zgadza się z oszacowaniem elementu wariancji.

Zauważ, że chociaż losowa procedura parowania jest intuicyjna i dydaktycznie przydatna, metoda zilustrowana przez @ Wolfganga jest w rzeczywistości o wiele mądrzejsza. W przypadku zestawu danych takiego jak ten o rozmiarze 100 * 100 liczba unikalnych par wewnątrzgrupowych (bez parowania własnego) wynosi 505 000 - duża, ale nie astronomiczna liczba - więc możemy całkowicie obliczyć korelację kompletnie wyczerpanego zestawu wszystkich możliwych par, zamiast konieczności losowego próbkowania z zestawu danych. Oto funkcja pobierania wszystkich możliwych par dla ogólnego przypadku z grupami dowolnej wielkości:

get_all_pairs <- function(df){
  # do this for every group and combine the results into a matrix
  do.call(rbind, by(df, df$group, function(group_df){
    # get all possible pairs of indices
    i <- expand.grid(seq(nrow(group_df)), seq(nrow(group_df)))
    # remove self-pairings
    i <- i[i[,1] != i[,2],]
    # return a 2-column matrix of the corresponding y-values
    cbind(group_df[i[,1], "y"], group_df[i[,2], "y"])
  }))
}

Teraz, jeśli zastosujemy tę funkcję do zestawu danych 100 * 100 i obliczymy korelację, otrzymamy:

cor(get_all_pairs(dat))
#           [,1]      [,2]
# [1,] 1.0000000 0.7504817
# [2,] 0.7504817 1.0000000

To, co dobrze zgadza się z pozostałymi dwoma oszacowaniami i w porównaniu z procedurą losowego parowania, jest znacznie szybsze do obliczenia, a także powinno być bardziej wydajnym oszacowaniem w sensie mniejszej wariancji.

Jake Westfall
źródło