Bakterien, die nach mehreren Oberflächenkontakten an den Fingern aufgenommen wurden: nicht normale Daten, wiederholte Messungen, gekreuzte Teilnehmer

9

Intro

Ich habe Teilnehmer, die unter zwei Bedingungen wiederholt kontaminierte Oberflächen mit E. coli berühren ( A = Handschuhe tragen, B = keine Handschuhe). Ich möchte wissen, ob es einen Unterschied zwischen der Menge an Bakterien auf ihren Fingerspitzen mit und ohne Handschuhe gibt, aber auch zwischen der Anzahl der Kontakte. Beide Faktoren sind innerhalb des Teilnehmers.

Experimentelle Methode:

Die Teilnehmer (n = 35) berühren jedes Quadrat einmal mit demselben Finger für maximal 8 Kontakte (siehe Abbildung a). a) Fingerkontakte mit 8 Oberflächen, b) KBE an den Fingern nach jedem Oberflächenkontakt

Ich wische dann den Finger des Teilnehmers ab und messe die Bakterien an der Fingerspitze nach jedem Kontakt. Sie berühren dann mit einem neuen Finger eine andere Anzahl von Oberflächen usw. von 1 bis 8 Kontakten (siehe Abbildung b).

Hier sind die realen Daten: reale Daten

Die Daten sind nicht normal, siehe unten unter Randverteilung von Bakterien | NumberContacts. x = Bakterien. Jede Facette hat eine andere Anzahl von Kontakten.

Geben Sie hier die Bildbeschreibung ein

MODELL

Versuch von lme4 :: glmer basierend auf Amöbenvorschlägen unter Verwendung von Gamma (link = "log") und Polynom für NumberContacts:

cfug<-glmer(CFU ~ Gloves + poly(NumberContacts,2) + (-1+NumberContacts|Participant),
            data=(K,CFU<4E5),
           family=Gamma(link="log")
            )
plot(cfug)

NB. Gamma (link = "inverse") wird nicht ausgeführt und sagt, dass die PIRLS-Halbierung die Abweichung nicht verringern konnte.

Ergebnisse:

Angepasst gegen Residuen für cfug Geben Sie hier die Bildbeschreibung ein

qqp (Resid (cfug))

Geben Sie hier die Bildbeschreibung ein

Frage:

Ist mein Glmer-Modell richtig definiert, um die zufälligen Effekte jedes Teilnehmers und die Tatsache zu berücksichtigen, dass jeder Experiment A gefolgt von Experiment B durchführt ?

Zusatz:

Autokorrelation scheint zwischen den Teilnehmern zu bestehen. Dies liegt wahrscheinlich daran, dass sie nicht am selben Tag getestet wurden und der Bakterienkolben mit der Zeit wächst und abnimmt. Ist das wichtig?

acf (CFU, lag = 35) zeigt eine signifikante Korrelation zwischen einem Teilnehmer und dem nächsten.

Geben Sie hier die Bildbeschreibung ein

HCAI
quelle
1
Sie können NumberContactseinen numerischen Faktor verwenden und quadratische / kubische Polynomterme einschließen. Oder schauen Sie sich Generalized Additive Mixed Models an.
Amöbe
1
@amoeba Danke für deine Hilfe. Alle Teilnehmer machten B (nicht behandelt), gefolgt von A (behandelt). Denken Sie, dass es andere grundlegende Probleme bei der Analyse gibt? Wenn ja, bin ich offen für weitere Antworten.
HCAI
1
Wenn ja, dann können Sie zufällige Wirkung des Handschuhs einschließen. Ich verstehe auch nicht, warum Sie zufällige Abschnitte entfernen und warum Sie nicht das gesamte Polynom 2. Grades in den zufälligen Teil aufnehmen. Und Sie können Handschuh * num Interaktion haben. Warum also nicht CFU ~ Gloves * poly(NumberContacts,2) + (Gloves * poly(NumberContacts,2) | Participant)oder so?
Amöbe
1
Oh, ich verstehe den Abschnitt, aber dann müssten Sie auch den festen Abschnitt unterdrücken. Auch für Nullkontakte sollten Sie keine KBE haben, aber mit Log-Link ist dies nicht sinnvoll. Und Sie haben bei 1 Kontakt nicht annähernd null KBE. Also würde ich das Abfangen nicht unterdrücken. Nicht konvergieren ist nicht gut, versuchen Sie, die Interaktion aus dem zufälligen Teil zu entfernen: CFU ~ Gloves * poly(NumberContacts,2) + (Gloves + poly(NumberContacts,2) | Participant)oder vielleicht die Handschuhe von dort zu entfernen CFU ~ Gloves * poly(NumberContacts,2) + (poly(NumberContacts,2) | Participant)...
Amöbe
1
Ich denke, es Gloves * poly(NumberContacts,2) + (poly(NumberContacts,2) | Participant)ist ein ziemlich anständiges Modell.
Amöbe

Antworten:

6

Einige Diagramme zum Erkunden der Daten

Unten sind acht, eine für jede Anzahl von Oberflächenkontakten, xy-Diagramme, die Handschuhe gegen keine Handschuhe zeigen.

Jedes Individuum ist mit einem Punkt versehen. Der Mittelwert sowie die Varianz und Kovarianz sind mit einem roten Punkt und der Ellipse angegeben (Mahalanobis-Abstand entspricht 97,5% der Bevölkerung).

14

Die kleine Korrelation zeigt, dass es tatsächlich einen zufälligen Effekt von den Individuen gibt (wenn es keinen Effekt von der Person gab, sollte es keine Korrelation zwischen den gepaarten Handschuhen und keinen Handschuhen geben). Dies ist jedoch nur ein kleiner Effekt, und eine Person kann unterschiedliche zufällige Effekte für "Handschuhe" und "keine Handschuhe" haben (z. B. kann die Person für alle unterschiedlichen Kontaktpunkte durchweg höhere / niedrigere Werte für "Handschuhe" als "keine Handschuhe" haben). .

xy Parzellen mit und ohne Handschuhe

Unterhalb des Diagramms befinden sich separate Diagramme für jede der 35 Personen. Die Idee dieses Diagramms ist es, zu sehen, ob das Verhalten homogen ist und welche Art von Funktion geeignet erscheint.

Beachten Sie, dass das "ohne Handschuhe" rot ist. In den meisten Fällen ist die rote Linie höher, mehr Bakterien für die Fälle "ohne Handschuhe".

Ich glaube, dass eine lineare Darstellung ausreichen sollte, um die Trends hier zu erfassen. Der Nachteil des quadratischen Diagramms besteht darin, dass die Koeffizienten schwieriger zu interpretieren sind (Sie werden nicht direkt sehen, ob die Steigung positiv oder negativ ist, da sowohl der lineare Term als auch der quadratische Term einen Einfluss darauf haben).

Aber was noch wichtiger ist, Sie sehen, dass die Trends zwischen den verschiedenen Individuen sehr unterschiedlich sind und es daher nützlich sein kann, einen zufälligen Effekt nicht nur für den Achsenabschnitt, sondern auch für die Steigung des Individuums hinzuzufügen.

Grundstücke für jeden Einzelnen

Modell

Mit dem Modell unten

  • Jedes Individuum erhält eine eigene angepasste Kurve (zufällige Effekte für lineare Koeffizienten).
  • yN(log(μ),σ2)log(y)N(μ,σ2)
  • Gewichte werden angewendet, weil die Daten heteroskedastisch sind. Die Variation ist zu den höheren Zahlen hin enger. Dies liegt wahrscheinlich daran, dass die Bakterienzahl eine gewisse Obergrenze hat und die Variation hauptsächlich auf eine fehlerhafte Übertragung von der Oberfläche auf den Finger zurückzuführen ist (= im Zusammenhang mit niedrigeren Zählungen). Siehe auch in den 35 Plots. Es gibt hauptsächlich einige wenige Personen, bei denen die Abweichung viel höher ist als bei den anderen. (Wir sehen auch größere Schwänze, Überdispersion, in den qq-Plots)
  • Es wird kein Intercept-Term verwendet und ein "Kontrast" -Term hinzugefügt. Dies geschieht, um die Interpretation der Koeffizienten zu erleichtern.

.

K    <- read.csv("~/Downloads/K.txt", sep="")
data <- K[K$Surface == 'P',]
Contactsnumber   <- data$NumberContacts
Contactscontrast <- data$NumberContacts * (1-2*(data$Gloves == 'U'))
data <- cbind(data, Contactsnumber, Contactscontrast)
m    <- lmer(log10CFU ~ 0 + Gloves + Contactsnumber + Contactscontrast + 
                        (0 + Gloves + Contactsnumber + Contactscontrast|Participant) ,
             data=data, weights = data$log10CFU)

Das gibt

> summary(m)
Linear mixed model fit by REML ['lmerMod']
Formula: log10CFU ~ 0 + Gloves + Contactsnumber + Contactscontrast + (0 +  
    Gloves + Contactsnumber + Contactscontrast | Participant)
   Data: data
Weights: data$log10CFU

REML criterion at convergence: 180.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0972 -0.5141  0.0500  0.5448  5.1193 

Random effects:
 Groups      Name             Variance  Std.Dev. Corr             
 Participant GlovesG          0.1242953 0.35256                   
             GlovesU          0.0542441 0.23290   0.03            
             Contactsnumber   0.0007191 0.02682  -0.60 -0.13      
             Contactscontrast 0.0009701 0.03115  -0.70  0.49  0.51
 Residual                     0.2496486 0.49965                   
Number of obs: 560, groups:  Participant, 35

Fixed effects:
                  Estimate Std. Error t value
GlovesG           4.203829   0.067646   62.14
GlovesU           4.363972   0.050226   86.89
Contactsnumber    0.043916   0.006308    6.96
Contactscontrast -0.007464   0.006854   -1.09

qqplot

Residuen

Code, um Diagramme zu erhalten

Chemometrie :: drawMahal-Funktion

# editted from chemometrics::drawMahal
drawelipse <- function (x, center, covariance, quantile = c(0.975, 0.75, 0.5, 
                                              0.25), m = 1000, lwdcrit = 1, ...) 
{
  me <- center
  covm <- covariance
  cov.svd <- svd(covm, nv = 0)
  r <- cov.svd[["u"]] %*% diag(sqrt(cov.svd[["d"]]))
  alphamd <- sqrt(qchisq(quantile, 2))
  lalpha <- length(alphamd)
  for (j in 1:lalpha) {
    e1md <- cos(c(0:m)/m * 2 * pi) * alphamd[j]
    e2md <- sin(c(0:m)/m * 2 * pi) * alphamd[j]
    emd <- cbind(e1md, e2md)
    ttmd <- t(r %*% t(emd)) + rep(1, m + 1) %o% me
#    if (j == 1) {
#      xmax <- max(c(x[, 1], ttmd[, 1]))
#      xmin <- min(c(x[, 1], ttmd[, 1]))
#      ymax <- max(c(x[, 2], ttmd[, 2]))
#      ymin <- min(c(x[, 2], ttmd[, 2]))
#      plot(x, xlim = c(xmin, xmax), ylim = c(ymin, ymax), 
#           ...)
#    }
  }
  sdx <- sd(x[, 1])
  sdy <- sd(x[, 2])
  for (j in 2:lalpha) {
    e1md <- cos(c(0:m)/m * 2 * pi) * alphamd[j]
    e2md <- sin(c(0:m)/m * 2 * pi) * alphamd[j]
    emd <- cbind(e1md, e2md)
    ttmd <- t(r %*% t(emd)) + rep(1, m + 1) %o% me
#    lines(ttmd[, 1], ttmd[, 2], type = "l", col = 2)
    lines(ttmd[, 1], ttmd[, 2], type = "l", col = 1, lty=2)  #
  }
  j <- 1
  e1md <- cos(c(0:m)/m * 2 * pi) * alphamd[j]
  e2md <- sin(c(0:m)/m * 2 * pi) * alphamd[j]
  emd <- cbind(e1md, e2md)
  ttmd <- t(r %*% t(emd)) + rep(1, m + 1) %o% me
#  lines(ttmd[, 1], ttmd[, 2], type = "l", col = 1, lwd = lwdcrit)
  invisible()
}

5 x 7 Grundstück

#### getting data
K <- read.csv("~/Downloads/K.txt", sep="")

### plotting 35 individuals

par(mar=c(2.6,2.6,2.1,1.1))
layout(matrix(1:35,5))

for (i in 1:35) {
  # selecting data with gloves for i-th participant
  sel <- c(1:624)[(K$Participant==i) & (K$Surface == 'P') & (K$Gloves == 'G')]
      # plot data
  plot(K$NumberContacts[sel],log(K$CFU,10)[sel], col=1,
       xlab="",ylab="",ylim=c(3,6))
      # model and plot fit
  m <- lm(log(K$CFU[sel],10) ~ K$NumberContacts[sel])
  lines(K$NumberContacts[sel],predict(m), col=1)

  # selecting data without gloves for i-th participant 
  sel <- c(1:624)[(K$Participant==i) & (K$Surface == 'P') & (K$Gloves == 'U')]
     # plot data 
  points(K$NumberContacts[sel],log(K$CFU,10)[sel], col=2)
     # model and plot fit
  m <- lm(log(K$CFU[sel],10) ~ K$NumberContacts[sel])
  lines(K$NumberContacts[sel],predict(m), col=2)
  title(paste0("participant ",i))
}

2 x 4 Grundstück

#### plotting 8 treatments (number of contacts)

par(mar=c(5.1,4.1,4.1,2.1))
layout(matrix(1:8,2,byrow=1))

for (i in c(1:8)) {
  # plot canvas
  plot(c(3,6),c(3,6), xlim = c(3,6), ylim = c(3,6), type="l", lty=2, xlab='gloves', ylab='no gloves')

  # select points and plot
  sel1 <- c(1:624)[(K$NumberContacts==i) & (K$Surface == 'P') & (K$Gloves == 'G')]
  sel2 <- c(1:624)[(K$NumberContacts==i) & (K$Surface == 'P') & (K$Gloves == 'U')]
  points(K$log10CFU[sel1],K$log10CFU[sel2])

  title(paste0("contact ",i))

  # plot mean
  points(mean(K$log10CFU[sel1]),mean(K$log10CFU[sel2]),pch=21,col=1,bg=2)

  # plot elipse for mahalanobis distance
  dd <- cbind(K$log10CFU[sel1],K$log10CFU[sel2])
  drawelipse(dd,center=apply(dd,2,mean),
            covariance=cov(dd),
            quantile=0.975,col="blue",
            xlim = c(3,6), ylim = c(3,6), type="l", lty=2, xlab='gloves', ylab='no gloves')
}
Sextus Empiricus
quelle
Vielen Dank Martijn, du hast die Dinge so klar erklärt. Tolle! Da das Kopfgeld abgelaufen ist, bevor ich es zuweisen konnte, möchte ich Ihnen sehr gerne einen separaten Betrag anbieten (ich werde mir jetzt ansehen, wie das geht). Ich habe jedoch einige Fragen: Erstens scheint die Transformation der Daten Denkschulen zu haben: Einige stimmen zu, andere stimmen vehement nicht zu. Warum ist es hier in Ordnung? Zweitens, warum erleichtert das Entfernen des zufälligen Abschnitts die Interpretation der Koeffizienten?
HCAI
(2) Ich denke, dass die Transformation in Ordnung ist, wenn man argumentieren könnte, dass es einen Prozess gibt, der die Transformation logisch macht (in der Tat widerstrebend zu transformieren, weil dadurch die Ergebnisse schön aussehen, kann dies als Datenmanipulation und falsche Darstellung der Ergebnisse sowie als Nichterhalten des Basiswerts angesehen werden Modell)
Sextus Empiricus
Ich sehe @Martijn, zumindest in der Biologie ist die Transformation durch log10 bei Bakterien üblich. Ich bin froh, das Kopfgeld zu geben, du hast es verdient. Würde es Ihnen etwas ausmachen, etwas näher darauf einzugehen, warum Sie diesen "Kontrastbegriff" verwenden?
HCAI
1
Bezüglich des Kontrasts Siehe hier stats.stackexchange.com/a/308644/164061 Sie haben die Freiheit, den Intercept-Term zu verschieben. Ein möglicherweise nützlicher Weg besteht darin, den Achsenabschnitt zwischen den beiden Kategorien festzulegen und den Effekt als Differenz zwischen den beiden Effekten (einer ist negativ, der andere positiv) relativ zu diesem mittleren Achsenabschnitt zu betrachten. (nicht, dass ich dafür eine Variable hinzufügen musste)
Sextus Empiricus
1
Idealerweise würden Sie die Behandlungen zufällig über die Zeit verteilen, so dass sich mögliche Auswirkungen aufgrund von Zeitschwankungen abschwächen würden. Aber ich sehe eigentlich nicht so viel Autokorrelation. Meinen Sie solche Sprünge wie bei Teilnehmer 5 zwischen 5 und 6 Kontakten, nach denen die Leitung wieder stabil ist? Ich denke, dass diese nicht so schlecht sind und höchstens Rauschen hinzufügen, aber Ihre Methode nicht stören (außer Signal / Rauschen niedrig zu machen). Sie können sicherer sein, wenn Sie im Laufe der Zeit keine systematische Änderung feststellen. Wenn Sie die Teilnehmer der Reihe nach verarbeitet haben, können Sie ihre mittlere KBE über die Zeit zeichnen.
Sextus Empiricus
2

Zur Frage , ob Gebrauch MASS:glmmPQLoder lme4:glmerfür Ihr Modell, ist mein Verständnis , dass diese beiden Funktionen werden das gleiche Modell passen (so lange , wie Sie die Modellgleichung, die Verteilung und die Verknüpfungsfunktion gleich eingestellt) , aber sie verwenden verschiedene Schätzmethoden die Passform zu finden. Ich könnte mich irren, aber mein Verständnis aus der Dokumentation ist, dass glmmPQLdie bestrafte Quasi-Wahrscheinlichkeit verwendet wird, wie in Wolfinger und O'Connell (1993) beschrieben , während glmerdie Gauß-Hermite-Quadratur verwendet wird. Wenn Sie sich darüber Sorgen machen, können Sie Ihr Modell mit beiden Methoden anpassen und überprüfen, ob sie dieselben Koeffizientenschätzungen liefern. Auf diese Weise können Sie sicherer sein, dass der Anpassungsalgorithmus zu den tatsächlichen MLEs der Koeffizienten konvergiert hat.


Sollte NumberContactsein kategorischer Faktor sein?

Diese Variable hat eine natürliche Reihenfolge, die aus Ihren Plots hervorgeht und eine reibungslose Beziehung zur Antwortvariablen aufweist, sodass Sie sie vernünftigerweise als numerische Variable behandeln können. Wenn Sie einbeziehen, werden factor(NumberContacts)Sie seine Form nicht einschränken und Sie werden nicht viele Freiheitsgrade verlieren. Sie können die Interaktion sogar nutzen, Gloves*factor(NumberContacts)ohne zu viele Freiheitsgrade zu verlieren. Es ist jedoch zu überlegen, ob die Verwendung einer Faktorvariablen eine Überanpassung der Daten beinhalten würde. Angesichts der Tatsache, dass Ihr Diagramm eine ziemlich glatte Beziehung aufweist, würde eine einfache lineare Funktion oder ein Quadrat gute Ergebnisse erzielen, ohne dass eine Überanpassung erforderlich ist.


Wie macht man Participanteine zufällige Steigung, aber keine Intercept-Variable?

Sie haben Ihre Antwortvariable bereits mithilfe einer logarithmischen Verknüpfungsfunktion auf eine Protokollskala gesetzt, sodass ein Abfangeffekt für Participanteinen multiplikativen Effekt auf die Antwort ergibt. Wenn Sie diesem eine zufällige Steigung geben NumberContactswürden, mit der er interagiert , hätte dies einen leistungsbasierten Effekt auf die Reaktion. Wenn Sie dies möchten, können Sie es erhalten, mit (~ -1 + NumberContacts|Participant)dem der Achsenabschnitt entfernt wird, aber eine Steigung basierend auf der Anzahl der Kontakte hinzugefügt wird.


Soll ich Box-Cox zum Transformieren meiner Daten verwenden? (zB Lambda = 0,779)

λ


Sollte ich Gewichte für die Varianz einschließen?

Sehen Sie sich zunächst Ihr Restdiagramm an, um festzustellen, ob Hinweise auf Heteroskedastizität vorliegen. Aufgrund der Diagramme, die Sie bereits aufgenommen haben, scheint es mir kein Problem zu sein, sodass Sie keine Gewichte für die Varianz hinzufügen müssen. Im Zweifelsfall können Sie Gewichte mit einer einfachen linearen Funktion hinzufügen und dann einen statistischen Test durchführen, um festzustellen, ob die Steigung der Gewichtung flach ist. Dies würde einen formalen Test der Heteroskedastizität bedeuten, der Ihnen ein Backup für Ihre Wahl geben würde.


Sollte ich Autokorrelation in einbeziehen NumberContacts?

Wenn Sie bereits einen zufälligen Effektbegriff für den Teilnehmer angegeben haben, ist es wahrscheinlich eine schlechte Idee, einen Autokorrelationsbegriff für die Anzahl der Kontakte hinzuzufügen. In Ihrem Experiment wird ein anderer Finger für eine unterschiedliche Anzahl von Kontakten verwendet, sodass Sie keine Autokorrelation für den Fall erwarten würden, in dem Sie den Teilnehmer bereits berücksichtigt haben. Das Hinzufügen eines Autokorrelationsbegriffs zusätzlich zum Teilnehmereffekt würde bedeuten, dass Sie glauben, dass eine bedingte Abhängigkeit zwischen dem Ergebnis verschiedener Finger besteht, basierend auf der Anzahl der Kontakte, selbst für einen bestimmten Teilnehmer.


Ben - Monica wieder einsetzen
quelle
Danke, das ist eine erstaunliche Antwort! Am Ende habe ich Gamma ausprobiert (link = "log") und der Glmer konvergiert ohne Beschwerde, Hurra! glmer (CFU ~ Handschuhe + Poly (NumberContacts, 2) + (-1 + NumberContacts | Teilnehmer), data = na.omit (Teilmenge (K, CFU <4.5e5 & Surface == "P")), family = Gamma ( link = "log")). QQplot ist meiner Meinung nach in Ordnung (nichts außerhalb von CIs), aber angepasst gegen Rediduals werden übertragen (siehe hinzugefügtes Bild, das hinzugefügt wurde, nachdem dieser Kommentar veröffentlicht wurde, falls er nicht übereinstimmt). Sollte ich mich zu sehr darum kümmern?
HCAI
1
QQ-Plot sieht für mich gut aus. Denken Sie auch daran, dass in einem GLM die Pearson-Residuen nicht unbedingt einer Normalverteilung folgen. Sieht so aus, als hätten Sie eine gute Analyse.
Ben - Reinstate Monica
1

In der Tat ist es vernünftig zu argumentieren, dass Messungen eines Teilnehmers nicht unabhängig von denen eines anderen Teilnehmers sind. Zum Beispiel neigen manche Menschen dazu, ihren Finger mit mehr (oder weniger) Kraft zu drücken, was sich auf alle ihre Messungen über jede Anzahl von Kontakten auswirken würde.

Daher wäre die 2-Wege-ANOVA mit wiederholten Messungen in diesem Fall ein akzeptables Modell.

Alternativ könnte man auch ein Modell mit gemischten Effekten participantals Zufallsfaktor anwenden . Dies ist eine fortschrittlichere und ausgefeiltere Lösung.

Mihael
quelle
Danke Mihael, du hast absolut Recht mit dem Druck. Hmm, ich habe hier über das Mixed-Effects-Modell rcompanion.org/handbook/I_09.html gelesen , bin mir aber nicht sicher über die Wechselwirkungen und verschachtelten Faktoren. Sind meine Faktoren verschachtelt?
HCAI
Ich möchte auch darauf hinweisen, dass die Daten nicht für jeden Kontakt normal verteilt sind . Sehen Sie sich daher die PQL-Modellierung (Penalized Quasi-Likelihood) an: ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/… . Denken Sie, dass dies eine gute Wahl ist?
HCAI