Bacterias detectadas en los dedos después de múltiples contactos en la superficie: datos no normales, medidas repetidas, participantes cruzados

9

Introducción

Tengo participantes que están tocando repetidamente superficies contaminadas con E. coli en dos condiciones ( A = usando guantes, B = sin guantes). Quiero saber si hay una diferencia entre la cantidad de bacterias en la punta de los dedos con y sin guantes, pero también entre la cantidad de contactos. Ambos factores están dentro del participante.

Método experimental:

Los participantes (n = 35) tocan cada cuadrado una vez con el mismo dedo para un máximo de 8 contactos (ver figura a). a) contactos de los dedos con 8 superficies, b) CFU en los dedos después de cada contacto de superficie

Luego froto el dedo del participante y mido las bacterias en la punta del dedo después de cada contacto. Luego usan un dedo nuevo para tocar un número diferente de superficies y así sucesivamente de 1 a 8 contactos (ver figura b).

Aquí están los datos reales : datos reales

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

ingrese la descripción de la imagen aquí

MODELO

Intentando desde lme4 :: glmer basado en las sugerencias de amebas 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)

NÓTESE BIEN. Gamma (enlace = "inverso") no se ejecutará diciendo que la reducción a la mitad PIRLS no pudo reducir la desviación.

Resultados:

Equipado vs residuales para cfug ingrese la descripción de la imagen aquí

qqp (resid (cfug))

ingrese la descripción de la imagen aquí

Pregunta:

¿Está mi modelo glmer correctamente definido para incorporar los efectos aleatorios de cada participante y el hecho de que todos hacen el experimento A seguido del experimento B ?

Adición:

La autocorrelación parece existir entre los participantes. Esto probablemente se deba a que no se analizaron el mismo día y el matraz de bacterias crece y disminuye con el tiempo. ¿Importa?

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

ingrese la descripción de la imagen aquí

HCAI
fuente
1
Puede usarlo NumberContactscomo factor numérico e incluir términos polinómicos cuadráticos / cúbicos. O busque en modelos mixtos aditivos generalizados.
ameba
1
@amoeba Gracias por tu ayuda. Todos los participantes hicieron B (sin guantes) seguido de A (con guantes). ¿Crees que hay otros problemas fundamentales con el análisis? Si es así, estoy abierto a cualquier otra respuesta.
HCAI
1
Si es así, puede incluir el efecto aleatorio del guante. Además, no entiendo por qué elimina la intercepción aleatoria y por qué no incluye todo el polinomio de segundo grado en la parte aleatoria. Y puedes tener interacción guante * num. Entonces, ¿por qué no CFU ~ Gloves * poly(NumberContacts,2) + (Gloves * poly(NumberContacts,2) | Participant)o algo así?
ameba
1
Oh, entiendo sobre la intercepción, pero entonces también necesitarías suprimir la intercepción fija. Además, para cero contactos debe tener cero CFU, pero con log-link esto no tiene sentido. Y no tiene nada cerca de cero UFC en 1 contacto. Entonces no suprimiría la intercepción. No converger no es bueno, intente eliminar la interacción de la parte aleatoria: CFU ~ Gloves * poly(NumberContacts,2) + (Gloves + poly(NumberContacts,2) | Participant)o tal vez quitar los guantes de allí CFU ~ Gloves * poly(NumberContacts,2) + (poly(NumberContacts,2) | Participant)...
ameba
1
Creo que Gloves * poly(NumberContacts,2) + (poly(NumberContacts,2) | Participant)es un modelo bastante decente.
ameba

Respuestas:

6

Algunas parcelas para explorar los datos.

A continuación hay ocho, uno para cada número de contactos de superficie, gráficos xy que muestran guantes versus no guantes.

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

14

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

xy parcelas de con y sin guantes

Debajo de la parcela hay parcelas separadas para cada una de las 35 personas. La idea de este gráfico es ver si el comportamiento es homogéneo y también qué tipo de función parece adecuada.

Tenga en cuenta 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 una trama lineal debería ser suficiente para capturar las tendencias aquí. La desventaja de la gráfica cuadrática es que los coeficientes serán más difíciles de interpretar (no verá directamente si la pendiente es positiva o negativa porque el término lineal y el término cuadrático influyen en esto).

Pero lo más importante es que ve que las tendencias difieren mucho entre los diferentes individuos y, por lo tanto, puede ser útil agregar un efecto aleatorio no solo para la intercepción, sino también la pendiente del individuo.

parcelas para cada individuo

Modelo

Con el modelo de abajo

  • Cada individuo obtendrá su propia curva ajustada (efectos aleatorios para coeficientes lineales).
  • yN(log(μ),σ2)log(y)N(μ,σ2)
  • Los pesos se aplican porque los datos son heteroscedasticos. La variación es más estrecha hacia los números más altos. Esto probablemente se deba a que el recuento de bacterias tiene cierto límite y la variación se debe principalmente a una falla en la transmisión de la superficie al dedo (= relacionado con recuentos más bajos). Ver también en las 35 parcelas. Hay principalmente unos pocos individuos para los cuales la variación es mucho mayor que los demás. (vemos también colas más grandes, sobredispersión, en las parcelas qq)
  • No se utiliza ningún término de intercepción y se agrega 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

qqplot

derechos residuales de autor

código para obtener parcelas

quimiometría :: función 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()
}

5 x 7 parcela

#### 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 parcela

#### 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')
}
Sexto empírico
fuente
Muchas gracias Martijn, has explicado las cosas muy claramente. ¡Asombroso! Como la recompensa terminó antes de que pudiera asignarla, me gustaría ofrecerle una cantidad separada (veré cómo hacerlo ahora). Sin embargo, tengo algunas dudas: en primer lugar, la transformación de los datos parece tener escuelas de pensamiento: algunos están de acuerdo y otros están totalmente en desacuerdo. ¿Por qué está bien aquí? En segundo lugar, ¿por qué eliminar la intersección aleatoria hace que los coeficientes sean más fáciles de interpretar?
HCAI
(2) Supongo que la transformación está bien cuando se podría argumentar que hay un proceso que hace que la transformación sea lógica (de hecho, la transformación a regañadientes porque hace que los resultados se vean bien puede verse como manipulación de datos y tergiversación de resultados, así como no obtener el subyacente modelo)
Sextus Empiricus
Veo @Martijn, al menos en biología la transformación por log10 es común para las bacterias. Estoy feliz de dar la recompensa, te lo mereces. ¿Te importaría elaborar un poco sobre por qué usas este "término de contraste" por favor?
HCAI
1
Con respecto al contraste Vea aquí stats.stackexchange.com/a/308644/164061 Usted tiene la libertad de mover el término de intercepción. Una forma posiblemente útil es establecer la intercepción entre las dos categorías y dejar que el efecto sea la diferencia entre los dos efectos (uno será negativo y el otro positivo) en relación con ese término medio de intercepción. (no es que tuviera que agregar una variable para esto)
Sextus Empiricus
1
Lo ideal sería que los tratamientos se distribuyan aleatoriamente a lo largo del tiempo, de modo que cualquier posible efecto debido a variaciones en el tiempo se nivele. Pero en realidad no veo tanta autocorrelación. ¿Te refieres a saltos como en el número 5 de contactos entre 5 y 6 participantes después de los cuales la línea vuelve a ser estable? Creo que estos no son tan malos y a lo sumo agregan ruido, pero no interfieren con su método (excepto hacer que la señal / ruido sea bajo). Puede estar más seguro cuando no ve un cambio sistemático en el tiempo. Si procesó a los participantes en orden, podría trazar su UFC media con el tiempo.
Sextus Empiricus
2

En cuanto a si usar MASS:glmmPQLo lme4:glmerpara su modelo, entiendo que ambas funciones se ajustarán al mismo modelo (siempre que configure la ecuación del modelo, la distribución y la función de enlace de la misma), pero utilizan diferentes métodos de estimación para encontrar el ajuste. Podría estar equivocado, pero mi comprensión de la documentación es que glmmPQLusa cuasi-verosimilitud penalizada como se describe en Wolfinger y O'Connell (1993) , mientras que glmerusa la cuadratura de Gauss-Hermite. Si le preocupa, puede ajustar su modelo con ambos métodos y verificar que brinden las mismas estimaciones de coeficientes y de esa manera tendrá una mayor confianza en que el algoritmo de ajuste ha convergido a los verdaderos MLE de los coeficientes.


¿Debería NumberContactsser un factor categórico?

Esta variable tiene un orden natural que parece tener una relación fluida con la variable de respuesta de sus gráficos, por lo que podría tratarla razonablemente como una variable numérica. Si fuera a incluir factor(NumberContacts), no limitará su forma y no perderá muchos grados de libertad. Incluso podría usar la interacción Gloves*factor(NumberContacts)sin perder demasiados grados de libertad. Sin embargo, vale la pena considerar si usar una variable de factor implicaría un ajuste excesivo de los datos. Dado que hay una relación bastante suave en su gráfica, una función lineal simple o cuadrática obtendría buenos resultados sin un ajuste excesivo.


¿Cómo se hace Participantuna pendiente aleatoria pero no se intercepta variable?

Ya ha puesto su variable de respuesta en una escala logarítmica utilizando una función de enlace logarítmico, por lo que un efecto de intercepción para Participantestá dando un efecto multiplicativo en la respuesta. Si le diera a esto una pendiente aleatoria con la que interactuaría NumberContacts, tendría un efecto basado en el poder en la respuesta. Si desea esto, puede obtenerlo con (~ -1 + NumberContacts|Participant)lo que eliminará la intersección pero agregará una pendiente en función del número de contactos.


¿Debo usar Box-Cox para transformar mis datos? (por ejemplo, lambda = 0.779)

λ


¿Debo incluir pesos para la varianza?

Comience mirando su parcela residual para ver si hay evidencia de heterocedasticidad. Según los gráficos que ya ha incluido, me parece que esto no es un problema, por lo que no necesita agregar ningún peso para la variación. En caso de duda, puede agregar pesos utilizando una función lineal simple y luego realizar una prueba estadística para ver si la pendiente de la ponderación es plana. Eso equivaldría a una prueba formal de heterocedasticidad, que le daría un respaldo para 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 agregar un término de autocorrelación en la cantidad 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 contabilizado al participante. Agregar un término de autocorrelación además del efecto participante significaría que cree que existe una dependencia condicional entre el resultado de diferentes dedos, en función del número de contactos, incluso para un participante dado.


Ben - Restablece a Monica
fuente
¡Gracias, esa es una respuesta increíble! Al final probé Gamma (link = "log") y el glmer converge sin quejas, ¡hurra! glmer (CFU ~ Gloves + poly (NumberContacts, 2) + (-1 + NumberContacts | Participant), data = na.omit (subconjunto (K, CFU <4.5e5 & Surface == "P")), family = Gamma ( link = "log")). QQplot, creo que está bien (nada fuera de los CI), pero los encajonados vs los rediduales están canalizando (vea la foto agregada después de publicar este comentario en caso de que no coincida). ¿Debería preocuparme demasiado por eso?
HCAI
1
QQ plot me parece bien. Además, recuerde que en un GLM los residuos de Pearson no siguen necesariamente una distribución normal. Parece que tienes un buen análisis.
Ben - Restablece a Mónica el
1

De hecho, es razonable argumentar que las mediciones tomadas de un participante no son independientes de las tomadas de otro participante. Por ejemplo, algunas personas pueden tender a presionar su dedo con más (o menos) fuerza, lo que afectaría todas sus mediciones en cada número de contactos.

Por lo tanto, el ANOVA de medidas repetidas de 2 vías sería un modelo aceptable para aplicar en este caso.

Alternativamente, uno también podría aplicar un modelo de efectos mixtos con participantun factor aleatorio. Esta es una solución más avanzada y más sofisticada.

Mihael
fuente
Gracias Mihael, tienes toda la razón sobre la presión. Hmm, estaba leyendo sobre el modelo de efectos mixtos aquí rcompanion.org/handbook/I_09.html pero no estoy seguro sobre las interacciones y los factores anidados. ¿Están anidados mis factores?
HCAI
También debo señalar que los datos no se distribuyen normalmente para cada contacto, así que he analizado el modelado de Cuasi-verosimilitud penalizada (PQL): ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/… . ¿Crees que esta es una buena opción?
HCAI