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

26

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.

Chris
quelle

Antworten:

21

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.

Henrik
quelle
12

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.

Mark Miller
quelle
10

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.). Bildbeschreibung hier eingeben

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)

Bildbeschreibung hier eingeben

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

alexforrence
quelle
1
+1. Natürlich ist die Varianz des Varianzschätzers umso größer, je kleiner die Anzahl der Subjekte ist. Aber ich denke nicht, dass es hier darauf ankommt. Die Frage ist, wie viele Themen ermöglichen es, vernünftige Ergebnisse zu erzielen? Wenn wir "unsinniges" Ergebnis als Null-Varianz-Ergebnis definieren, passiert dies in Ihrer Simulation ziemlich oft mit n = 5 oder weniger. Ausgehend von n = 6 oder n = 7 erhalten Sie fast nie eine exakte 0-Varianzschätzung, dh das Modell konvergiert zu einer nicht entarteten Lösung. Meine Schlussfolgerung wäre, dass n = 6 grenzwertig akzeptabel ist.
Amöbe sagt Reinstate Monica
1
Übrigens steht dies im Einklang mit rpubs.com/bbolker/4187 .
Amöbe sagt Reinstate Monica
8

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.

StasK
quelle
1
Ich gehe davon aus, dass die Antwort auf die Frage gewissermaßen lautet: "Wie lang ist ein Stück Schnur?". Aber ich würde nicht viel Vertrauen in die Schätzung eines Mittelwerts oder einer Varianz aus einer Stichprobe von weniger als 15-20 setzen, also würde dieselbe Faustregel nicht für Stufen eines zufälligen Effekts gelten. Ich habe noch nie gesehen, dass jemand in Längsschnittstudien die Themen-ID als festen und zufälligen Effekt einbezog. Ist dies eine gängige Praxis?
Chris
Abgesehen von einer kleinen Anzahl von Probanden in einem gemischten Modell werden ihre zufälligen Effekte nicht beobachtet, sodass Sie sie aus den Daten herausfiltern müssen, und Sie benötigen wahrscheinlich relativ mehr Daten, um dies verlässlich zu tun, als bei der bloßen Schätzung des Mittelwerts und der Varianz, wenn alles beobachtet wird. Also 42 gegen 15-20 :). Ich denke, ich meinte zufällige Steigungen, da Sie in den Subjekt-IDs, die nur als zufällige Effekte dargestellt werden, richtig sind, sonst werden sie nicht identifiziert. Ökonomen glauben im Übrigen nicht an zufällige Effekte und veröffentlichen fast ausschließlich das, was sie als "feste Effekte" bezeichnen, dh subjektinterne Schätzungen.
StasK
2
+1 @StasK für eine sehr gute Antwort auf eine Frage, die sehr schwer zu beantworten ist. Ich denke, es gibt einen Anflug von unnötigem Sarkasmus, und Sie könnten in Erwägung ziehen, Ihre Antwort so zu bearbeiten, dass Sie dem OP ein wenig mehr Respekt entgegenbringen.
Michael R. Chernick
@Michael, Sie haben wahrscheinlich Recht, dass dies eine launische Antwort ist, und wahrscheinlich unnötigerweise. Das OP akzeptierte jedoch die Antwort, die sie hören wollten, und erhielt daher eine Entschließung zu diesem Thema. Eine ernstere Antwort würde entweder auf einen guten Simulationsnachweis oder auf eine asymptotische Analyse höherer Ordnung hinweisen, aber mir sind solche Referenzen leider nicht bekannt.
StasK
3
Für das, was es wert ist, denke ich, geht es bei der magischen Zahl "42" nicht darum, ob zufällige Effekte gerechtfertigt sind, sondern ob man davonkommen kann, ohne sich um endlich große Korrekturen zu sorgen (z. B. Gedanken über effektive Nenner-Freiheitsgrade / Kenward-Roger-Korrekturen / andere ähnliche Ansätze).
Ben Bolker
7

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.

Tom Wenseleers
quelle
4

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.

Charles
quelle
Wurde aus einem Grund, der mit dem technischen Inhalt zu tun hat, dafür gestimmt?
N Brouwer
Mit welcher Art von Daten arbeiten Sie? Ich bin mir nicht sicher, warum Sie überrascht sind, dass das Modell mit 12-16 Messungen pro Person konvergiert. Ich kann die Verzerrung in den resultierenden Modellen nicht kommentieren, aber ich habe nie Probleme mit der Konvergenz in lme4gemischten Modellen gehabt und sie häufig mit ähnlichen Stichprobengrößen wie das OP ausgeführt (ich arbeite auch mit Biologie-Datensätzen).
RTbecard
1

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.

m_e_s
quelle