Bactéries ramassées sur les doigts après plusieurs contacts de surface: données non normales, mesures répétées, participants croisés

Intro

J'ai des participants qui touchent à plusieurs reprises des surfaces contaminées avec E. coli dans deux conditions ( A = porter des gants, B = pas de gants). Je veux savoir s'il y a une différence entre la quantité de bactéries au bout de leurs doigts avec et sans gants mais aussi entre le nombre de contacts. Les deux facteurs sont intra-participants.

Méthode expérimentale:

Les participants (n = 35) touchent chaque carré une fois avec le même doigt pour un maximum de 8 contacts (voir figure a).

Je tamponne ensuite le doigt du participant et mesure les bactéries sur le bout du doigt après chaque contact. Ils utilisent ensuite un nouveau doigt pour toucher un nombre différent de surfaces et ainsi de suite de 1 à 8 contacts (voir figure b).

Voici les vraies données: les vraies données

Les données ne sont pas normales, voir la distribution marginale des bactéries | NumberContacts ci-dessous. x = bactérie. Chaque facette est un nombre différent de contacts.

MODÈLE

Essayer de lme4 :: glmer basé sur les suggestions d'amibe en utilisant Gamma (link = "log") et polynôme pour NumberContacts:

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

NB. Le gamma (link = "inverse") ne fonctionnera pas en disant que la réduction de moitié des pas PIRLS n'a pas réussi à réduire la déviance.

Résultats:

Ajusté vs résiduels pour cfug

qqp (resid (cfug))

Question:

Mon modèle glmer est-il correctement défini pour incorporer les effets aléatoires de chaque participant et le fait que tout le monde fait à la fois l'expérience A puis l'expérience B ?

Une addition:

L'autocorrélation semble exister entre les participants. C'est probablement parce qu'ils n'ont pas été testés le même jour et que le flacon de bactéries se développe et diminue au fil du temps. Est-ce que ça importe?

acf (CFU, lag = 35) montre une corrélation significative entre un participant et le suivant.

Quelques parcelles pour explorer les données

Ci-dessous, huit, un pour chaque nombre de contacts de surface, des tracés xy montrant des gants par rapport à l'absence de gants.

Chaque individu est tracé avec un point. La moyenne et la variance et la covariance sont indiquées par un point rouge et l'ellipse (distance de Mahalanobis correspondant à 97,5% de la population).

Vous pouvez voir que les effets ne sont que faibles par rapport à la répartition de la population. La moyenne est plus élevée pour «pas de gants» et la moyenne change un peu plus haut pour plus de contacts de surface (ce qui peut être significatif). Mais l'effet n'est que de petite taille ( réduction globale du log ), et il y a beaucoup d'individus pour qui il y a en fait un plus grand nombre de bactéries avec les gants.$\frac{1}{4}$

La petite corrélation montre qu'il y a en effet un effet aléatoire des individus (s'il n'y a pas eu d'effet de la personne alors il ne devrait pas y avoir de corrélation entre les gants appariés et pas de gants). Mais ce n'est qu'un petit effet et un individu peut avoir des effets aléatoires différents pour les «gants» et «pas de gants» (par exemple, pour tous les points de contact différents, l'individu peut avoir systématiquement un nombre plus élevé / plus faible de «gants» que de «pas de gants») .

Le graphique ci-dessous présente des graphiques distincts pour chacun des 35 individus. L'idée de cette intrigue est de voir si le comportement est homogène et aussi de voir quel type de fonction semble convenir.

Notez que le «sans gants» est en rouge. Dans la plupart des cas, la ligne rouge est plus élevée, plus de bactéries pour les cas «sans gants».

Je crois qu'un tracé linéaire devrait être suffisant pour saisir les tendances ici. L'inconvénient du tracé quadratique est que les coefficients seront plus difficiles à interpréter (vous ne verrez pas directement si la pente est positive ou négative car le terme linéaire et le terme quadratique ont tous deux une influence sur cela).

Mais plus important encore, vous voyez que les tendances diffèrent beaucoup entre les différents individus et qu'il peut donc être utile d'ajouter un effet aléatoire non seulement pour l'interception, mais aussi pour la pente de l'individu.

Modèle

Avec le modèle ci-dessous

• Chaque individu obtiendra sa propre courbe ajustée (effets aléatoires pour les coefficients linéaires).
• Le modèle utilise des données transformées en logarithme et correspond à un modèle linéaire régulier (gaussien). Dans les commentaires, l'amibe a mentionné qu'un lien de journal n'est pas lié à une distribution lognormale. Mais c'est différent. est différent delog ( y ) ∼ N ( μ , σ 2$y \sim N(\log(\mu),\sigma^2)$$\log(y) \sim N(\mu,\sigma^2)$
• Des pondérations sont appliquées car les données sont hétéroscédastiques. La variation est plus étroite vers les nombres plus élevés. Cela est probablement dû au fait que le nombre de bactéries a un certain plafond et que la variation est principalement due à une transmission défaillante de la surface au doigt (= liée à un nombre inférieur). Voir aussi dans les 35 graphiques. Il y a principalement quelques individus pour lesquels la variation est beaucoup plus élevée que les autres. (on voit aussi des queues plus grosses, une surdispersion, dans les parcelles qq)
• Aucun terme d'interception n'est utilisé et un terme de «contraste» est ajouté. Ceci est fait pour rendre les coefficients plus faciles à interpréter.

.

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)

Cela donne

> 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

code pour obtenir des tracés

chimiométrie :: fonction 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 parcelle

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

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

Tout d'abord, bon travail sur votre graphique; il donne une représentation claire des données, vous pouvez donc déjà voir le type de motif dans les données en fonction du nombre de contacts et de l'utilisation ou de l'absence de gants. En regardant ce graphique, je pense que vous obtiendrez de bons résultats avec un modèle log-polynomial de base, avec des effets aléatoires pour les participants. Le modèle que vous avez choisi semble raisonnable, mais vous pouvez également envisager d'ajouter un terme quadratique pour le nombre de contacts.$^\dagger$

Quant à savoir s'il faut utiliser MASS:glmmPQLou lme4:glmerpour votre modèle, je crois comprendre que ces deux fonctions s'adapteront au même modèle (tant que vous définissez l'équation, la distribution et la fonction de lien du modèle de la même manière), mais elles utilisent différentes méthodes d'estimation pour trouver l'ajustement. Je peux me tromper, mais d'après ce que je comprends de la documentation, on glmmPQLutilise la quasi-vraisemblance pénalisée comme décrit dans Wolfinger et O'Connell (1993) , alors que l'on glmerutilise la quadrature de Gauss-Hermite. Si cela vous inquiète, vous pouvez ajuster votre modèle avec les deux méthodes et vérifier qu'elles donnent les mêmes estimations de coefficient et de cette façon, vous aurez une plus grande confiance que l'algorithme d'ajustement a convergé vers les vrais MLE des coefficients.

Devrait NumberContactsêtre un facteur catégorique?

Cette variable a un ordre naturel qui apparaît à partir de vos graphiques pour avoir une relation fluide avec la variable de réponse, vous pouvez donc raisonnablement la traiter comme une variable numérique. Si vous deviez l'inclure, factor(NumberContacts)vous ne contraindriez pas sa forme et vous ne perdriez pas beaucoup de degrés de liberté. Vous pouvez même utiliser l'interaction Gloves*factor(NumberContacts)sans perdre trop de degrés de liberté. Cependant, il convient de se demander si l'utilisation d'une variable factorielle impliquerait un surajustement des données. Étant donné qu'il existe une relation assez fluide dans votre graphique, une simple fonction linéaire ou quadratique obtiendrait de bons résultats sans surajustement.

Comment faire Participantune pente aléatoire mais pas une variable d'interception?

Vous avez déjà placé votre variable de réponse sur une échelle logarithmique en utilisant une fonction de lien logarithmique, donc un effet d'interception pour Participantdonne un effet multiplicatif sur la réponse. Si vous deviez lui donner une pente aléatoire en interaction, NumberContactscela aurait un effet basé sur la puissance sur la réponse. Si vous le souhaitez, vous pouvez l'obtenir avec (~ -1 + NumberContacts|Participant)ce qui supprimera l'interception mais ajoutera une pente en fonction du nombre de contacts.

Dois-je utiliser Box-Cox pour transformer mes données? (par exemple lambda = 0,779)

En cas de doute, essayez d'ajuster un modèle avec cette transformation et voyez comment elle se compare à d'autres modèles en utilisant des statistiques de qualité d'ajustement appropriées. Si vous allez utiliser cette transformation, il serait préférable de laisser le paramètre tant que paramètre libre et de le laisser être estimé dans le cadre de votre modèle, plutôt que de pré-spécifier une valeur.$\lambda$

Dois-je inclure des pondérations pour la variance?

Commencez par regarder votre tracé résiduel pour voir s'il existe des preuves d'hétéroscédasticité. Sur la base des graphiques que vous avez déjà inclus, il me semble que ce n'est pas un problème, vous n'avez donc pas besoin d'ajouter de pondération pour la variance. En cas de doute, vous pouvez ajouter des poids à l'aide d'une simple fonction linéaire, puis effectuer un test statistique pour voir si la pente de la pondération est plate. Cela équivaudrait à un test formel d'hétéroscédasticité, ce qui vous donnerait une sauvegarde de votre choix.

Dois-je inclure l'autocorrélation dans NumberContacts?

Si vous avez déjà inclus un terme à effet aléatoire pour le participant, ce serait probablement une mauvaise idée d'ajouter un terme d'auto-corrélation sur le nombre de contacts. Votre expérience utilise un doigt différent pour différents nombres de contacts, vous ne vous attendez donc pas à une autocorrélation dans le cas où vous avez déjà comptabilisé le participant. L'ajout d'un terme d'autocorrélation en plus de l'effet participant signifierait que vous pensez qu'il existe une dépendance conditionnelle entre le résultat de différents doigts, en fonction du nombre de contacts, même pour un participant donné.

$^\dagger$ Votre graphique montre bien les relations, mais vous pouvez l'améliorer esthétiquement en ajoutant un titre et des sous-titres et en lui donnant de meilleures étiquettes d'axe. Vous pouvez également simplifier votre légende en supprimant son titre et en remplaçant «Oui» par «Gants» et «Non» par «Pas de gants».

En effet, il est raisonnable de soutenir que les mesures prises par un participant ne sont pas indépendantes de celles prises par un autre participant. Par exemple, certaines personnes pourraient avoir tendance à appuyer sur leur doigt avec plus (ou moins) de force, ce qui affecterait toutes leurs mesures sur chaque nombre de contacts.

L'ANOVA à mesures répétées bidirectionnelles serait donc un modèle acceptable à appliquer dans ce cas.

Alternativement, on pourrait également appliquer un modèle à effets mixtes avec participantcomme facteur aléatoire. Il s'agit d'une solution plus avancée et plus sophistiquée.