Bakterie zbierane na palcach po wielu kontaktach powierzchniowych: dane nienormalne, powtarzane pomiary, skrzyżowane osoby

Wprowadzenie

Mam uczestników, którzy wielokrotnie dotykają skażonych powierzchni E. coli w dwóch warunkach ( A = w rękawiczkach, B = bez rękawiczek). Chcę wiedzieć, czy istnieje różnica między ilością bakterii na opuszkach palców w rękawicach i bez rękawiczek, ale także między liczbą kontaktów. Oba czynniki są wewnątrz uczestników.

Metoda eksperymentalna:

Uczestnicy (n = 35) dotykają każdego kwadratu tym samym palcem maksymalnie przez 8 kontaktów (patrz rysunek a).

Następnie przecieram palec uczestnika i po każdym kontakcie mierzę bakterie na opuszku palca. Następnie używają nowego palca do dotykania innej liczby powierzchni i od 1 do 8 kontaktów (patrz rysunek b).

Oto prawdziwe dane: prawdziwe dane

Dane są nienormalne, więc zobacz marginalny rozkład bakterii | NumberContacts poniżej. x = bakterie. Każdy aspekt to inna liczba kontaktów.

MODEL

Próbowanie z lme4 :: glmer na podstawie sugestii ameby za pomocą Gamma (link = "log") i wielomian dla 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”) nie uruchomi się, mówiąc, że skok o połowę PIRLS nie zmniejszył dewiacji.

Wyniki:

Dopasowane vs resztki dla cfug

qqp (resid (cfug))

Pytanie:

Czy mój model blasku jest odpowiednio zdefiniowany, aby uwzględnić losowe efekty każdego uczestnika oraz fakt, że każdy wykonuje zarówno eksperyment A, jak i eksperyment B ?

Dodanie:

Wydaje się, że istnieje autokorelacja między uczestnikami. Jest tak prawdopodobnie dlatego, że nie zostały przetestowane tego samego dnia, a kolba bakterii rośnie i maleje z czasem. Czy to ma znaczenie?

acf (CFU, lag = 35) pokazuje istotną korelację między jednym uczestnikiem a drugim.

Niektóre wykresy do eksploracji danych

Poniżej znajduje się osiem, po jednym dla każdej liczby styków powierzchniowych, wykresy xy pokazujące rękawice w porównaniu z brakiem rękawic.

Każda osoba jest wykreślona kropką. Średnia, wariancja i kowariancja są oznaczone czerwoną kropką i elipsą (odległość Mahalanobisa odpowiadająca 97,5% populacji).

Widać, że efekty są niewielkie w porównaniu do rozprzestrzeniania się populacji. Średnia jest wyższa dla „bez rękawiczek”, a średnia zmienia się nieco wyżej w przypadku większej liczby kontaktów na powierzchni (co można wykazać jako znaczące). Ale efekt jest tylko trochę w rozmiarze (ogólnie obniżenie log), i istnieje wiele osób, bo kto tam jest w rzeczywistości wyższe bakterie liczyć z rękawic.$\frac{1}{4}$

Mała korelacja pokazuje, że rzeczywiście istnieje efekt losowy od poszczególnych osób (jeśli nie wystąpiłby efekt od tej osoby, wówczas nie powinno być żadnej korelacji między sparowanymi rękawiczkami i bez rękawiczek). Jest to jednak tylko niewielki efekt, a osoba może mieć różne losowe skutki dla „rękawiczek” i „bez rękawiczek” (np. Dla wszystkich różnych punktów kontaktowych dana osoba może mieć stale wyższą / niższą liczbę dla „rękawiczek” niż „bez rękawiczek”) .

Poniżej wykresu są osobne wykresy dla każdej z 35 osobników. Ideą tego wykresu jest sprawdzenie, czy zachowanie jest jednorodne, a także sprawdzenie, jaki rodzaj funkcji wydaje się odpowiedni.

Pamiętaj, że „bez rękawiczek” ma kolor czerwony. W większości przypadków czerwona linia jest wyższa, więcej bakterii w przypadkach „bez rękawiczek”.

Uważam, że liniowa fabuła powinna wystarczyć do uchwycenia trendów tutaj. Wadą wykresu kwadratowego jest to, że współczynniki będą trudniejsze do interpretacji (nie zobaczysz bezpośrednio, czy nachylenie jest dodatnie czy ujemne, ponieważ wpływ na to mają zarówno element liniowy, jak i kwadratowy).

Ale co ważniejsze, widać, że trendy różnią się znacznie między poszczególnymi osobami, dlatego może być użyteczny dodanie losowego efektu nie tylko dla przechwytywania, ale także nachylenia jednostki.

Model

Z poniższym modelem

• Każda osoba otrzyma własną krzywą (efekty losowe dla współczynników liniowych).
• Model wykorzystuje dane przekształcone w log i pasuje do zwykłego (gaussowskiego) modelu liniowego. W komentarzach amoeba wspomniała, że ​​link do dziennika nie odnosi się do rozkładu logarytmicznego. Ale to jest inne.$y \sim N(\log(\mu),\sigma^2)$ jest różny od $\log(y) \sim N(\mu,\sigma^2)$
• Wagi są stosowane, ponieważ dane są heteroskedastyczne. Różnica jest węższa w stosunku do wyższych liczb. Wynika to prawdopodobnie z tego, że liczba bakterii ma pewien pułap, a zmienność wynika głównie z nieudanej transmisji z powierzchni na palec (= związane z mniejszą liczbą). Zobacz także na 35 działkach. Istnieje głównie kilka osobników, dla których zróżnicowanie jest znacznie wyższe niż u innych. (widzimy także większe ogony, nadmierną dyspersję, na wykresach qq)
• Nie stosuje się terminu przechwytującego i dodaje się termin „kontrastowy”. Ma to na celu ułatwienie interpretacji współczynników.

.

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)

To daje

> 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

kod do uzyskania wykresów

chemometrics :: funkcja drawMahal

# 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()
}

Działka 5 x 7

#### 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))
}

Działka 2 x 4

#### 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')
}

Po pierwsze, dobra praca na wykresie; daje wyraźną reprezentację danych, dzięki czemu można już zobaczyć rodzaj wzorca w danych na podstawie liczby kontaktów oraz użycia lub braku rękawiczek.$^\dagger$ Patrząc na ten wykres, myślę, że uzyskasz dobre wyniki dzięki prostemu modelowi log-wielomianu, z losowymi efektami dla uczestników. Wybrany model wygląda rozsądnie, ale możesz również rozważyć dodanie kwadratowego terminu dla liczby kontaktów.

Jeśli chodzi o to, czy użyć MASS:glmmPQLlub lme4:glmerdla twojego modelu, rozumiem, że obie te funkcje będą pasowały do ​​tego samego modelu (o ile ustawisz równanie modelu, rozkład i funkcję połączenia tak samo), ale używają różnych metod szacowania, aby znaleźć dopasowanie. Mogę się mylić, ale rozumiem z dokumentacji, że glmmPQLwykorzystuje karę quasi-prawdopodobieństwa, jak opisano w Wolfinger i O'Connell (1993) , podczas gdy glmerużywa kwadratury Gaussa-Hermite'a. Jeśli martwisz się o to, możesz dopasować swój model do obu metod i sprawdzić, czy dają one takie same oszacowania współczynników, w ten sposób zyskasz większą pewność, że algorytm dopasowywania zbliżył się do prawdziwych MLE współczynników.

Czy powinien to NumberContactsbyć czynnik kategoryczny?

Ta zmienna ma naturalny porządek, który pojawia się na twoich wykresach, aby mieć płynny związek ze zmienną odpowiedzi, więc możesz rozsądnie traktować ją jako zmienną numeryczną. Jeśli miałbyś to uwzględnić factor(NumberContacts), nie ograniczysz jego formy i nie stracisz wielu stopni swobody. Możesz nawet skorzystać z interakcji Gloves*factor(NumberContacts)bez utraty zbyt wielu stopni swobody. Warto jednak zastanowić się, czy użycie zmiennej czynnikowej wiązałoby się z nadmiernym dopasowaniem danych. Biorąc pod uwagę, że na wykresie istnieje dość gładka zależność, prosta funkcja liniowa lub kwadratowa uzyskałaby dobre wyniki bez nadmiernego dopasowania.

Jak utworzyć Participantlosowe nachylenie, ale nie przechwycić zmiennej?

Umieściłeś już swoją zmienną odpowiedzi w skali logarytmicznej za pomocą logarytmicznej funkcji link, więc efekt przechwytywania dla Participantdaje efekt multiplikatywny dla odpowiedzi. Gdyby nadać temu losowe nachylenie wchodzące w interakcje NumberContacts, miałoby to zależny od mocy wpływ na odpowiedź. Jeśli chcesz, możesz to uzyskać, dzięki (~ -1 + NumberContacts|Participant)czemu usuniesz przechwytywanie, ale doda nachylenie na podstawie liczby kontaktów.

Czy powinienem używać Box-Cox do transformacji moich danych? (np. lambda = 0,779)

W razie wątpliwości spróbuj dopasować model do tej transformacji i zobacz, jak wypada on w porównaniu z innymi modelami, korzystając z odpowiednich statystyk dobroci dopasowania. Jeśli zamierzasz użyć tej transformacji, lepiej pozostaw parametr$\lambda$ jako darmowy parametr i pozwól, aby był szacowany jako część modelu, zamiast wcześniejszego określania wartości.

Czy powinienem dołączyć wagi dla wariancji?

Zacznij od spojrzenia na resztkowy wykres, aby sprawdzić, czy istnieją dowody na heteroscedastyczność. Na podstawie fabuł, które już załączyłeś, wydaje mi się, że to nie jest problem, więc nie musisz dodawać żadnych wag dla wariancji. W razie wątpliwości można dodać ciężary za pomocą prostej funkcji liniowej, a następnie wykonać test statystyczny, aby sprawdzić, czy nachylenie ważenia jest płaskie. Oznaczałoby to formalny test heteroscedastyczności, który dałby ci pewne wsparcie dla twojego wyboru.

Czy powinienem włączyć autokorelację NumberContacts?

Jeśli podałeś już uczestnikowi pojęcie efektu losowego, prawdopodobnie złym pomysłem byłoby dodanie terminu autokorelacji do liczby kontaktów. W twoim eksperymencie używa się innego palca dla różnych liczb kontaktów, więc nie spodziewałbyś się autokorelacji w przypadku, gdy już rozliczałeś uczestnika. Dodanie terminu autokorelacji oprócz efektu uczestnika oznaczałoby, że uważasz, że istnieje zależność warunkowa między wynikami różnych palców, oparta na liczbie kontaktów, nawet dla danego uczestnika.

$^\dagger$Twój wykres dobrze pokazuje relacje, ale możesz go poprawić estetycznie, dodając informacje o tytule i podtytule i nadając mu lepsze etykiety osi. Możesz także uprościć swoją legendę, usuwając jej tytuł i zmieniając „Tak” na „Rękawiczki” i „Nie” na „Brak rękawiczek”.

Rzeczywiście uzasadnione jest twierdzenie, że pomiary wykonane od jednego uczestnika nie są niezależne od pomiarów wykonanych od innego uczestnika. Na przykład niektóre osoby mogą naciskać palcem z większą (lub mniejszą) siłą, co wpłynęłoby na wszystkie ich pomiary dla każdej liczby kontaktów.

Zatem dwukierunkowa ANOVA z powtarzanymi pomiarami byłaby akceptowalnym modelem do zastosowania w tym przypadku.

Alternatywnie można również zastosować model efektów mieszanych participantjako czynnik losowy. Jest to bardziej zaawansowane i bardziej wyrafinowane rozwiązanie.