@Wolfgang hat schon eine tolle Antwort gegeben. Ich möchte es ein wenig erweitern, um zu zeigen, dass Sie auch zu dem geschätzten ICC von 0,75 in seinem Beispieldatensatz gelangen können, indem Sie den intuitiven Algorithmus der zufälligen Auswahl vieler Paare von Werten buchstäblich implementieren - wobei die Mitglieder jedes Paares von der stammen gleiche Gruppe - und dann einfach ihre Korrelation berechnen. Und dann kann das gleiche Verfahren leicht auf Datensätze mit Gruppen beliebiger Größe angewendet werden, wie ich auch zeigen werde.y
Zuerst laden wir den Datensatz von @ Wolfgang (hier nicht gezeigt). Definieren wir nun eine einfache R-Funktion, die ein data.frame nimmt und ein einzelnes zufällig ausgewähltes Beobachtungspaar aus derselben Gruppe zurückgibt:
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"])
}
Hier ist ein Beispiel dafür, was wir erhalten, wenn wir diese Funktion 10 Mal auf dem Datensatz von @ Wolfgang aufrufen:
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
Um nun den ICC abzuschätzen, rufen wir diese Funktion einfach mehrmals auf und berechnen dann die Korrelation zwischen den beiden Spalten.
random_pairs <- replicate(100000, get_random_pair(dat))
cor(t(random_pairs))
# [,1] [,2]
# [1,] 1.0000000 0.7493072
# [2,] 0.7493072 1.0000000
Das gleiche Verfahren kann ohne Änderungen auf Datensätze mit Gruppen beliebiger Größe angewendet werden. Erstellen wir beispielsweise einen Datensatz, der aus 100 Gruppen mit jeweils 100 Beobachtungen besteht, wobei der wahre ICC wie im Beispiel von @ Wolfgang auf 0,75 festgelegt ist.
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")
Wenn wir den ICC anhand der Varianzkomponenten eines gemischten Modells schätzen, erhalten wir:
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
Und wenn wir das zufällige Paarungsverfahren anwenden, erhalten wir
random_pairs <- replicate(100000, get_random_pair(dat))
cor(t(random_pairs))
# [,1] [,2]
# [1,] 1.0000000 0.7503004
# [2,] 0.7503004 1.0000000
was eng mit der Varianzkomponentenschätzung übereinstimmt.
Beachten Sie, dass das Zufallspaarungsverfahren zwar intuitiv und didaktisch nützlich ist, die von @Wolfgang veranschaulichte Methode jedoch viel intelligenter ist. Für einen Datensatz wie diesen mit der Größe 100 * 100 beträgt die Anzahl der eindeutigen gruppeninternen Paarungen (ohne Selbstpaarungen) 505.000 - eine große, aber keine astronomische Zahl -, sodass wir die Korrelation durchaus berechnen können der vollständig erschöpften Menge aller möglichen Paarungen, anstatt zufällig aus dem Datensatz zu stichproben. Hier ist eine Funktion, um alle möglichen Paarungen für den allgemeinen Fall mit Gruppen beliebiger Größe abzurufen:
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"])
}))
}
Wenn wir nun diese Funktion auf den 100 * 100-Datensatz anwenden und die Korrelation berechnen, erhalten wir:
cor(get_all_pairs(dat))
# [,1] [,2]
# [1,] 1.0000000 0.7504817
# [2,] 0.7504817 1.0000000
Was gut mit den anderen beiden Schätzungen übereinstimmt und mit dem Zufallspaarungsverfahren verglichen wird, ist viel schneller zu berechnen und sollte auch eine effizientere Schätzung im Sinne einer geringeren Varianz sein.
