Bacterias recogidas en los dedos tras múltiples contactos con la superficie: datos no normales, medidas repetidas, participantes cruzados

Question

Introducción

Tengo participantes que tocan repetidamente superficies contaminadas con E. coli en dos condiciones ( A \=llevar guantes, B \=sin guantes). Quiero saber si hay diferencia entre la cantidad de bacterias en las yemas de los dedos con y sin guantes, pero también entre el número de contactos. Ambos factores son intraparticipantes.

Método experimental:

Los participantes (n=35) tocan cada cuadrado una vez con el mismo dedo hasta un máximo de 8 contactos (véase la figura a).

A continuación, froto el dedo del participante y mido las bacterias en la yema del dedo después de cada contacto. A continuación, utilizan un nuevo dedo para tocar un número diferente de superficies y así sucesivamente de 1 a 8 contactos (véase la figura b).

Estos son los datos reales: datos reales

Los datos no son normales, por lo que véase la distribución marginal de bacterias|NúmeroContactos a continuación. x=bacterias. Cada faceta es un número diferente de contactos.

MODELO

Tratando de lme4::glmer basado en las sugerencias de amoeba usando Gamma(link="log") y polinomio para NumberContacts:

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

NB. Gamma(link="inversa") no se ejecuta porque el análisis por pasos de PIRLS no consigue reducir la desviación.

Resultados:

Ajuste frente a residuos para cfug

qqp(resid(cfug))

Pregunta:

¿Está bien definido mi modelo glmer para incorporar los efectos aleatorios de cada participante y el hecho de que todos hagan ambos experimentos? A seguido de un experimento B ?

Adición:

Parece existir autocorrelación entre los participantes. Esto se debe probablemente a que no se hicieron las pruebas el mismo día y a que el frasco de bacterias crece y disminuye con el tiempo. ¿Importa?

acf(UFC,lag=35) muestra una correlación significativa entre un participante y el siguiente.

Answer 1

Algunos gráficos para explorar los datos

A continuación se muestran ocho gráficos x-y, uno para cada número de contactos con la superficie, que muestran los guantes frente a la ausencia de guantes.

Cada individuo se representa con un punto. La media, la varianza y la covarianza se indican con un punto rojo y la elipse (distancia de Mahalanobis correspondiente al 97,5% de la población).

Se puede ver que los efectos son pequeños en comparación con la extensión de la población. La media es más alta para "sin guantes" y la media cambia un poco más arriba para más contactos superficiales (lo que puede demostrarse que es significativo). Pero el efecto es muy pequeño (en general, a $\frac{1}{4}$ reducción logarítmica), y hay muchos individuos para los que en realidad hay un mayor recuento de bacterias con los guantes.

La pequeña correlación muestra que, efectivamente, existe un efecto aleatorio de los individuos (si no hubiera un efecto de la persona, entonces no debería haber correlación entre los guantes emparejados y sin guantes). Pero es sólo un efecto pequeño y un individuo puede tener diferentes efectos aleatorios para "guantes" y "sin guantes" (por ejemplo, para todos los diferentes puntos de contacto el individuo puede tener recuentos consistentemente más altos/bajos para "guantes" que para "sin guantes").

Debajo se muestran parcelas separadas para cada uno de los 35 individuos. La idea de este gráfico es ver si el comportamiento es homogéneo y también ver qué tipo de función parece adecuada.

Nótese que el "sin guantes" está en rojo. En la mayoría de los casos la línea roja es más alta, más bacterias para los casos 'sin guantes'.

Creo que un gráfico lineal debería bastar para captar las tendencias. La desventaja del gráfico cuadrático es que los coeficientes van a ser más difíciles de interpretar (no vas a ver directamente si la pendiente es positiva o negativa porque tanto el término lineal como el cuadrático influyen en ello).

Pero lo más importante es que se ve que las tendencias difieren mucho entre los distintos individuos y, por lo tanto, puede ser útil añadir un efecto aleatorio no sólo para el intercepto, sino también para la pendiente del individuo.

Modelo

Con el modelo siguiente

• Cada individuo tendrá su propia curva ajustada (efectos aleatorios para coeficientes lineales).
• El modelo utiliza datos transformados logarítmicamente y se ajusta con un modelo lineal regular (gaussiano). En los comentarios, ameba mencionó que un enlace logarítmico no se relaciona con una distribución lognormal. Pero esto es diferente. $y \sim N(\log(\mu),\sigma^2)$ es diferente de $\log(y) \sim N(\mu,\sigma^2)$
• Se aplican ponderaciones porque los datos son heteroscedásticos. La variación es más estrecha hacia las cifras más altas. Esto se debe probablemente a que el recuento de bacterias tiene un cierto techo y la variación se debe sobre todo a la falta de transmisión de la superficie al dedo (= relacionado con recuentos más bajos). Véase también el gráfico 35. Hay principalmente unos pocos individuos para los que la variación es mucho mayor que para los demás. (vemos también colas más grandes, sobredispersión, en los gráficos qq)
• No se utiliza ningún término de intercepción y se añade un término de "contraste". Esto se hace para que los coeficientes sean más fáciles de interpretar.

.

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)

Esto da

> 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

código para obtener parcelas

función chemometrics::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()
}

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

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

Answer 2

En primer lugar, buen trabajo con tu gráfico; ofrece una representación clara de los datos, por lo que ya se puede ver el tipo de patrón en los datos basado en el número de contactos y el uso o ausencia de guantes. $^\dagger$ Observando este gráfico, creo que obtendrá buenos resultados con un modelo log-polinomial básico, con efectos aleatorios para los participantes. El modelo que ha elegido parece razonable, pero también podría considerar añadir un término cuadrático para el número de contactos.

En cuanto a la utilización de MASS:glmmPQL o lme4:glmer para su modelo, tengo entendido que ambas funciones se ajustarán al mismo modelo (siempre que establezca la ecuación del modelo, la distribución y la función de enlace de la misma manera), pero utilizan diferentes métodos de estimación para encontrar el ajuste. Puede que me equivoque, pero según la documentación entiendo que glmmPQL utiliza la cuasi verosimilitud penalizada descrita en Wolfinger y O'Connell (1993) mientras que glmer utiliza la cuadratura Gauss-Hermite. Si le preocupa, puede ajustar el modelo con ambos métodos y comprobar que dan las mismas estimaciones de los coeficientes; de este modo tendrá más confianza en que el algoritmo de ajuste ha convergido a los verdaderos MLE de los coeficientes.

---

Debería NumberContacts ser un factor categórico?

Esta variable tiene un orden natural que, según los gráficos, parece tener una relación suave con la variable de respuesta, por lo que podría tratarse razonablemente como una variable numérica. Si incluyera factor(NumberContacts) entonces no limitarás su forma y no perderás muchos grados de libertad. Incluso podría utilizar la interacción Gloves*factor(NumberContacts) sin perder demasiados grados de libertad. Sin embargo, conviene plantearse si el uso de una variable factorial implicaría un ajuste excesivo de los datos. Dado que el gráfico muestra una relación bastante suave, una función lineal simple o cuadrática daría buenos resultados sin sobreajustar los datos.

---

¿Cómo se hace Participant ¿una variable aleatoria de pendiente pero no de intercepto?

Ya ha colocado su variable de respuesta en una escala logarítmica utilizando una función de enlace logarítmica, por lo que un efecto de intercepción para Participant está dando un efecto multiplicativo en la respuesta. Si se le diera una pendiente aleatoria interactuando con NumberContacts entonces tendría un efecto basado en el poder en la respuesta. Si quiere esto, puede conseguirlo con (~ -1 + NumberContacts|Participant) que eliminará el intercepto pero añadirá una pendiente basada en el número de contactos.

---

¿Debo utilizar Box-Cox para transformar mis datos? (por ejemplo, lambda=0,779)

En caso de duda, pruebe a ajustar un modelo con esta transformación y vea cómo se compara con otros modelos utilizando las estadísticas de bondad de ajuste adecuadas. Si va a utilizar esta transformación, sería mejor dejar el parámetro $\lambda$ como parámetro libre y deje que se estime como parte de su modelo, en lugar de especificar previamente un valor.

---

¿Debo incluir ponderaciones para la varianza?

Empiece por observar el gráfico de residuos para ver si hay indicios de heteroscedasticidad. Basándome en los gráficos que ya has incluido, me parece que no es un problema, por lo que no necesitas añadir ninguna ponderación para la varianza. En caso de duda, podría añadir ponderaciones utilizando una función lineal simple y, a continuación, realizar una prueba estadística para comprobar si la pendiente de la ponderación es plana. Eso equivaldría a una prueba formal de heteroscedasticidad, lo que le daría cierto respaldo a su elección.

---

¿Debo incluir la autocorrelación en NumberContacts ?

Si ya ha incluido un término de efecto aleatorio para el participante, probablemente sería una mala idea añadir un término de autocorrelación en el número de contactos. Su experimento utiliza un dedo diferente para diferentes números de contactos, por lo que no esperaría autocorrelación para el caso en el que ya ha tenido en cuenta al participante. Añadir un término de autocorrelación además del efecto del participante significaría que usted cree que existe una dependencia condicional entre el resultado de diferentes dedos, basada en el número de contactos, incluso para un participante determinado.

---

$^\dagger$ Tu gráfico es bueno para mostrar las relaciones, pero podrías mejorarlo estéticamente añadiendo un título e información de subtítulos y dándole mejores etiquetas para los ejes. También podrías simplificar tu leyenda eliminando su título y cambiando "Sí" por "Guantes" y "No" por "Sin guantes".

Answer 3

De hecho, es razonable argumentar que las mediciones tomadas de uno de los participantes no son independientes de las que se toman de otro participante. Por ejemplo, algunas personas podrían tiende a presionar su dedo con más (o menos) de la fuerza, que afectaría a todas sus mediciones a través de cada número de contactos.

Así que el 2 vías de medidas repetidas ANOVA sería aceptable modelo a aplicar en este caso.

Alternativamente, también se podría aplicar un modelo de efectos mixtos con participant como un factor aleatorio. Esta es una de las más avanzadas y una solución más sofisticada.