Was ist die empfohlene Mindestanzahl von Gruppen für einen Zufallseffektfaktor?

Ich verwende ein gemischtes Modell in R( lme4), um einige Messwiederholungsdaten zu analysieren. Ich habe eine Reaktionsvariable (Fasergehalt von Kot) und 3 feste Effekte (Körpermasse usw.). Meine Studie hat nur 6 Teilnehmer mit jeweils 16 Wiederholungsmessungen (obwohl zwei nur 12 Wiederholungen haben). Die Versuchspersonen sind Eidechsen, denen in verschiedenen "Behandlungen" verschiedene Futterkombinationen verabreicht wurden.

Meine Frage ist: Kann ich die Betreff-ID als zufälligen Effekt verwenden?

Ich weiß, dass dies die übliche Vorgehensweise bei Longitudinal-Mixed-Effects-Modellen ist, um der zufällig ausgewählten Natur der Probanden und der Tatsache Rechnung zu tragen, dass Beobachtungen innerhalb von Probanden enger korrelieren als diejenigen zwischen Probanden. Wenn Sie die Subjekt-ID jedoch als zufälligen Effekt behandeln, müssen Sie einen Mittelwert und eine Varianz für diese Variable schätzen.

• Reicht dies aus, um den Mittelwert und die Varianz genau zu charakterisieren, da ich nur 6 Probanden habe (6 Stufen dieses Faktors)?

• Hilft in dieser Hinsicht die Tatsache, dass ich für jedes Thema einige wiederholte Messungen habe (ich sehe nicht, wie wichtig das ist)?

• Wenn ich die Betreff-ID nicht als zufälligen Effekt verwenden kann, kann ich dann durch Einbeziehen der ID als festen Effekt überprüfen, ob ich wiederholte Messungen durchgeführt habe.

Bearbeiten: Ich möchte nur klarstellen, dass wenn ich sage "Kann ich" Betreff-ID als zufälligen Effekt verwenden, ich meine "ist es eine gute Idee zu". Ich weiß, dass ich das Modell mit einem Faktor von nur 2 Stufen ausstatten kann, aber das wäre doch nicht zu rechtfertigen. Ich frage, wann es sinnvoll wird, Themen als zufällige Effekte zu behandeln. In der Literatur wird anscheinend davon ausgegangen, dass 5-6 Stufen eine Untergrenze sind. Es scheint mir, dass die Schätzungen des Mittelwerts und der Varianz des Zufallseffekts nicht sehr genau wären, bis es 15+ Faktorstufen gäbe.

Kurze Antwort: Ja, Sie können ID als Zufallseffekt mit 6 Stufen verwenden.

Etwas längere Antwort: In der GLMM-FAQ von @ BenBolker steht unter der Überschrift " Soll ich Faktor xxx als fest oder zufällig behandeln? ":

Ein Punkt von besonderer Relevanz für die 'moderne' gemischte Modellschätzung (anstelle der 'klassischen' Methode der Momentschätzung) ist, dass es aus praktischen Gründen eine vernünftige Anzahl von Zufallseffektniveaus (z. B. Blöcke) geben muss - mehr als Mindestens 5 oder 6.

Sie befinden sich also an der unteren Grenze, aber auf der rechten Seite.

Um die Mindestanzahl von Gruppen für ein Mehrebenenmodell zu ermitteln, habe ich mir das Buch Datenanalyse mit Regression und Mehrebenen / Hierarchische Modelle von Gelman und Hill (2007) angesehen.

Sie scheinen sich mit diesem Thema in Kapitel 11, Abschnitt 5 (Seite 247) zu befassen, wo sie schreiben, dass bei <5 Gruppen Mehrebenenmodelle in der Regel weniger als klassische Modelle hinzufügen. Sie scheinen jedoch zu schreiben, dass die Anwendung eines mehrstufigen Modells nur ein geringes Risiko birgt.

Dieselben Autoren scheinen zu diesem Thema in Kapitel 12, Abschnitt 9 (Seiten 275-276) zurückzukehren. Dort schreiben sie, dass Ratschläge zur Mindestanzahl von Gruppen für ein Mehrebenenmodell falsch sind. Dort heißt es wiederum, dass Mehrebenenmodelle bei geringer Gruppenanzahl häufig nur wenig zu klassischen Modellen beitragen. Sie schreiben jedoch auch, dass Mehrebenenmodelle nicht schlechter sein sollten als No-Pooling-Regression (wobei No-Pooling zu bedeuten scheint, dass Gruppenindikatoren in der klassischen Regression verwendet werden).

Auf den Seiten 275-276 haben die Autoren einen speziellen Unterabschnitt für den Fall von einer oder zwei Gruppen (z. B. männlich gegen weiblich). Hier schreiben sie, dass sie das Modell typischerweise in klassischer Form ausdrücken. Sie geben jedoch an, dass die mehrstufige Modellierung auch mit nur einer oder zwei Gruppen nützlich sein kann. Sie schreiben, dass sich die Mehrebenenmodellierung mit einer oder zwei Gruppen auf die klassische Regression reduziert.

Mein Eindruck davon ist, dass die klassische Regression ein Ende eines Kontinuums von Modellen ist, dh ein Sonderfall eines Mehrebenenmodells.

Auf der Grundlage der obigen Ausführungen gehe ich davon aus, dass die klassische Regression und die Mehrebenenmodellierung bei nur zwei Gruppen nahezu identische Schätzungen liefern und die Verwendung von Mehrebenenmodellen mit nur einer, zwei, drei, vier, fünf oder sechs Gruppen in Ordnung ist.

Ich werde versuchen, diese Antwort in Zukunft mit RCode und einem kleinen Datensatz zu modifizieren , der Schätzungen vergleicht, die mit beiden Ansätzen bei Verwendung von zwei Gruppen erhalten wurden.

Für das, was es wert ist, habe ich eine kleine Simulationsstudie durchgeführt, um die Stabilität der Varianzschätzung für ein relativ einfaches LMM (unter Verwendung des sleepstudydurch verfügbaren Datensatzes lme4) zu untersuchen. Der erste Weg generiert alle möglichen Subjektkombinationen für die ngroupsAnzahl der Subjekte und passt das Modell für jede mögliche Kombination an. Die Sekunde nimmt einige gelegentliche Teilmengen von Themen.

library(lme4)
library(ggplot2)
library(tidyr)

m0 <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy,
           control = lmerControl(optimizer = "nloptwrap"))
# set the number of factor levels
ngroups <- 3:18 
# generate all possible combinations
combos <- lapply(X = ngroups, 
                 FUN = function(x) combn(unique(sleepstudy$Subject), x)) 

# allocate output (sorry, this code is entirely un-optimized)
out <- list(matrix(NA, ncol(combos[[1]]), 1), matrix(NA, ncol(combos[[2]]), 1),
            matrix(NA, ncol(combos[[3]]), 1), matrix(NA, ncol(combos[[4]]), 1),
            matrix(NA, ncol(combos[[5]]), 1), matrix(NA, ncol(combos[[6]]), 1),
            matrix(NA, ncol(combos[[7]]), 1), matrix(NA, ncol(combos[[8]]), 1),
            matrix(NA, ncol(combos[[9]]), 1), matrix(NA, ncol(combos[[10]]), 1),
            matrix(NA, ncol(combos[[11]]), 1), matrix(NA, ncol(combos[[12]]), 1),
            matrix(NA, ncol(combos[[13]]), 1), matrix(NA, ncol(combos[[14]]), 1),
            matrix(NA, ncol(combos[[15]]), 1), matrix(NA, ncol(combos[[16]]), 1))
# took ~ 2.5 hrs on my laptop, commented out for safety
#system.time(for(ii in 1:length(combos)) {
#    for(jj in 1:ncol(combos[[ii]])) {
#    sls <- sleepstudy[sleepstudy$Subject %in% combos[[ii]][,jj],]
#    out[[ii]][jj] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
#        }
#    })

# pad with zeros, not all were equal
# from http://stackoverflow.com/questions/11148429/r-convert-asymmetric-list-to-matrix-number-of-elements-in-each-sub-list-diffe
max.len <- max(sapply(out, length))
corrected.list <- lapply(out, function(x) {c(x, rep(NA, max.len - length(x)))})
mat <- do.call(rbind, corrected.list)
mat <- data.frame(t(mat))
names(mat) <- paste0('s',3:18)
mat <- gather(mat, run, value)

ggplot(mat, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

Die gepunktete schwarze Linie ist die ursprüngliche Punktschätzung der Varianz, und die Facetten repräsentieren unterschiedliche Anzahlen von Themen ( s3als Gruppen von drei Themen, s4als vier, usw.).

Und der alternative Weg:

ngroups <- 3:18
reps <- 500
out2<- matrix(NA, length(ngroups), reps)

for (ii in 1:length(ngroups)) {
    for(j in 1:reps) {
        sls <- sleepstudy[sleepstudy$Subject %in% sample(unique(sleepstudy$Subject), ngroups[i], replace = FALSE),]
        out2[i,j] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
    }
}
out2 <- data.frame(t(out2))
names(out2) <- paste0('s',3:18)
out2 <- gather(out2, run, value)

ggplot(out2, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

Es scheint (für dieses Beispiel jedenfalls), dass sich die Varianz nicht wirklich stabilisiert, bis es mindestens 14 Probanden gibt, wenn nicht später.

Angrist und Pischkes "Mostly Harmless Econometrics" haben einen Abschnitt mit dem Titel "Weniger als 42 Cluster", in dem sie halb im Scherz sagen:

Nach der Annahme, dass die Antwort auf das Leben, das Universum und alles 42 ist, lautet die Frage daher: Wie viele Cluster reichen aus, um mit der Standard-Clusteranpassung [ähnlich dem Varianzschätzer in GEE] zuverlässig zu schließen?

Mein Ausbilder für Ökonometrie beantwortete Fragen wie Ihre wie folgt: "Amerika ist ein freies Land. Sie können tun, was Sie wollen. Wenn Sie jedoch Ihre Arbeit veröffentlichen möchten, müssen Sie in der Lage sein, das zu verteidigen, was Sie getan haben. " Mit anderen Worten, Sie werden wahrscheinlich in der Lage sein, R- oder Stata- oder HLM- oder Mplus- oder SAS PROC GLIMMIX-Code mit 6 Betreffs auszuführen (und zu diesen alternativen Paketen zu wechseln, wenn eines Ihrer Wahl dies nicht ausführt) sehr schwierige Zeit, diesen Ansatz zu verteidigen und asymptotische Tests zu rechtfertigen.

Ich glaube, dass standardmäßig das Einschließen einer Variablen als zufällige Steigung auch das Einschließen dieses Effekts als festen Effekt impliziert, und dass Sie durch viele Syntaxrahmen springen müssen, wenn Sie dies nur als zufälligen Effekt mit dem Mittelwert von haben möchten Null. Das ist eine vernünftige Entscheidung, die die Softwareentwickler für Sie getroffen haben.

Sie können auch ein Bayes'sches Mischmodell verwenden. In diesem Fall wird die Unsicherheit bei der Schätzung der zufälligen Effekte bei der Berechnung der glaubwürdigen 95% -Vorhersageintervalle vollständig berücksichtigt. Das neue R-Paket brmsund die neue R- Funktion ermöglichen brmzum Beispiel einen sehr einfachen Übergang von einem lme4frequentistischen gemischten Modell zu einem Bayes-Modell, da die Syntax nahezu identisch ist.

Ich würde kein Zufallseffektmodell mit nur 6 Ebenen verwenden. Modelle, die einen 6-stufigen Zufallseffekt verwenden, können manchmal mit vielen statistischen Programmen ausgeführt werden und bieten manchmal unvoreingenommene Schätzungen, aber:

1. Ich denke, es gibt einen willkürlichen Konsens in der statistischen Gemeinschaft, dass 10-20 die Mindestanzahl ist. Wenn Sie möchten, dass Ihre Forschung veröffentlicht wird, sollten Sie nach einer Zeitschrift ohne statistische Überprüfung suchen (oder Ihre Entscheidung in einer recht raffinierten Sprache begründen).
2. Bei so wenigen Clustern ist es wahrscheinlich, dass die Varianz zwischen den Clustern schlecht geschätzt wird. Eine schlechte Schätzung der Varianz zwischen den Clustern führt normalerweise zu einer schlechten Schätzung des Standardfehlers der interessierenden Koeffizienten. (Modelle mit zufälligen Effekten hängen von der Anzahl der Cluster ab, die theoretisch unendlich werden).
3. Oft konvergieren die Modelle einfach nicht. Haben Sie versucht, Ihr Modell zu betreiben? Ich würde mit nur 12-16 Maßnahmen pro Thema überraschen, dass die Modelle zusammenlaufen. Als ich es geschafft habe, diese Art von Modell zu konvergieren, hatte ich Hunderte von Messungen pro Cluster.

Dieses Problem wird in den meisten Standardlehrbüchern auf diesem Gebiet behandelt, und Sie haben es in Ihrer Frage sozusagen angesprochen. Ich glaube nicht, dass ich Ihnen neue Informationen gebe.

Die ursprüngliche Frage ist lange her, aber ich dachte, ich könnte ein paar Punkte hinzufügen, die für die Modellauswahl relevant sind.

1 - Solange das Modell identifiziert ist (dh Sie haben Freiheitsgrade im Parameterraum), sollten Sie VERSUCHEN können, das Modell anzupassen. Je nach Optimierungsmethode kann das Modell konvergieren oder nicht. Auf keinen Fall würde ich versuchen, mehr als 1 oder 2 zufällige Effekte und definitiv nicht mehr als 1 ebenenübergreifende Interaktion einzubeziehen. Im konkreten Fall des hier vorgestellten Problems reicht die Gruppengröße 6 möglicherweise nicht aus, um hinreichend genaue Schätzungen vorzunehmen, wenn wir eine Wechselwirkung zwischen eidechsenspezifischen Merkmalen (z. B. Alter, Größe usw.) und Behandlungs- / Messmerkmalen vermuten.

2 - Wie in einigen Antworten erwähnt, kann Konvergenz ein Problem sein. Ich habe jedoch die Erfahrung gemacht, dass sozialwissenschaftliche Daten aufgrund von Messproblemen ein großes Konvergenzproblem aufweisen, während Biowissenschaften und insbesondere biochemische Wiederholungsmessungen wesentlich kleinere Standardfehler aufweisen. Es hängt alles vom Prozess der Datengenerierung ab. Bei sozialen und wirtschaftlichen Daten müssen wir auf verschiedenen Abstraktionsebenen arbeiten. In biologischen und chemischen und mit Sicherheit astronomischen Daten ist ein Messfehler weniger ein Problem.