Bactérias apanhadas nos dedos após vários contatos superficiais: dados não normais, medidas repetidas, participantes cruzados

Introdução

Tenho participantes que tocam repetidamente superfícies contaminadas com E. coli em duas condições ( A = usando luvas, B = sem luvas). Quero saber se há uma diferença entre a quantidade de bactérias na ponta dos dedos com e sem luvas, mas também entre o número de contatos. Ambos os fatores são participantes.

Método experimental:

Os participantes (n = 35) tocam em cada quadrado uma vez com o mesmo dedo para um máximo de 8 contatos (veja a figura a).

Em seguida, esfrego o dedo do participante e meço as bactérias na ponta do dedo após cada contato. Eles então usam um novo dedo para tocar em um número diferente de superfícies e assim por diante de 1 a 8 contatos (consulte a figura b).

Aqui estão os dados reais : dados reais

Os dados não são normais, portanto veja a distribuição marginal de bactérias | NumberContacts abaixo. x = bactérias. Cada faceta é um número diferente de contatos.

MODELO

Tentando em lme4 :: glmer com base nas sugestões da ameba usando Gamma (link = "log") e polinômio 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 = "inverso") não será executado dizendo que a etapa PIRLS não conseguiu reduzir o desvio.

Resultados:

Ajustado vs resíduos para cfug

qqp (resid (cfug))

Questão:

Meu modelo de glmer está definido adequadamente para incorporar os efeitos aleatórios de cada participante e o fato de que todo mundo faz o experimento A seguido pelo experimento B ?

Adição:

A autocorrelação parece existir entre os participantes. Provavelmente, porque não foram testados no mesmo dia e o frasco de bactérias cresce e diminui com o tempo. Isso importa?

acf (UFC, lag = 35) mostra uma correlação significativa entre um participante e o próximo.

Alguns gráficos para explorar os dados

Abaixo estão oito, um para cada número de contatos de superfície, gráficos xy mostrando luvas versus sem luvas.

Cada indivíduo é plotado com um ponto. A média e variância e covariância são indicadas com um ponto vermelho e a elipse (distância de Mahalanobis correspondente a 97,5% da população).

Você pode ver que os efeitos são apenas pequenos em comparação com a expansão da população. A média é mais alta para 'sem luvas' e a média muda um pouco mais para obter mais contatos de superfície (que podem ser significativos). Mas o efeito é apenas de tamanho pequeno ( redução geral de um ), e há muitos indivíduos para quem há realmente uma contagem mais alta de bactérias com as luvas.$\frac{1}{4}$

A pequena correlação mostra que há realmente um efeito aleatório dos indivíduos (se não houve um efeito da pessoa, não deve haver correlação entre as luvas emparelhadas e as sem luvas). Mas é apenas um efeito pequeno e um indivíduo pode ter efeitos aleatórios diferentes para 'luvas' e 'sem luvas' (por exemplo, para todos os diferentes pontos de contato, o indivíduo pode ter contagens sempre mais altas / baixas de 'luvas' do que 'sem luvas') .

A plotagem abaixo contém plotagens separadas para cada um dos 35 indivíduos. A idéia desse gráfico é verificar se o comportamento é homogêneo e também ver que tipo de função parece adequada.

Observe que o 'sem luvas' está em vermelho. Na maioria dos casos, a linha vermelha é mais alta, mais bactérias para os casos 'sem luvas'.

Eu acredito que um gráfico linear deve ser suficiente para capturar as tendências aqui. A desvantagem do gráfico quadrático é que os coeficientes serão mais difíceis de interpretar (você não verá diretamente se a inclinação é positiva ou negativa porque o termo linear e o termo quadrático influenciam isso).

Mais importante, porém, você vê que as tendências diferem muito entre os diferentes indivíduos e, portanto, pode ser útil adicionar um efeito aleatório não apenas para a interceptação, mas também para a inclinação do indivíduo.

Modelo

Com o modelo abaixo

• Cada indivíduo terá sua própria curva ajustada (efeitos aleatórios para coeficientes lineares).
• O modelo usa dados transformados em log e se ajusta a um modelo linear regular (gaussiano). Nos comentários, a ameba mencionou que um link de log não está relacionado a uma distribuição lognormal. Mas isso é diferente. é diferente delog ( y ) ∼ N ( μ , σ 2$y \sim N(\log(\mu),\sigma^2)$$\log(y) \sim N(\mu,\sigma^2)$
• Os pesos são aplicados porque os dados são heterocedásticos. A variação é mais estreita em relação aos números mais altos. Provavelmente, isso ocorre porque a contagem de bactérias tem algum limite e a variação se deve principalmente à falha na transmissão da superfície para o dedo (= relacionada a contagens mais baixas). Veja também nas 35 parcelas. Existem principalmente alguns indivíduos para os quais a variação é muito maior que os outros. (vemos também caudas maiores, superdispersão, nos gráficos de qq)
• Nenhum termo de interceptação é usado e um termo de 'contraste' é adicionado. Isso é feito para facilitar a interpretação dos coeficientes.

.

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)

Isto dá

> 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 obter parcelas

quimiometria :: função 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()
}

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

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

$^\dagger$

Quanto a usar MASS:glmmPQLou lme4:glmerpara o seu modelo, meu entendimento é que essas duas funções se encaixam no mesmo modelo (desde que você defina a equação do modelo, a distribuição e a função de link da mesma forma), mas usam métodos de estimativa diferentes para encontrar o ajuste. Eu poderia estar enganado, mas meu entendimento na documentação é que glmmPQLusa quase-probabilidade penalizada, conforme descrito em Wolfinger e O'Connell (1993) , enquanto glmerusa quadratura de Gauss-Hermite. Se você estiver preocupado com isso, poderá ajustar seu modelo com ambos os métodos e verificar se eles fornecem as mesmas estimativas de coeficiente e, dessa forma, você terá maior confiança de que o algoritmo de ajuste convergiu para os verdadeiros MLEs dos coeficientes.

Deve NumberContactsser um fator categórico?

Essa variável possui uma ordem natural que aparece nos seus gráficos para ter um relacionamento suave com a variável de resposta, para que você possa tratá-la razoavelmente como uma variável numérica. Se você incluir factor(NumberContacts), não restringirá sua forma e não perderá muitos graus de liberdade. Você pode até usar a interação Gloves*factor(NumberContacts)sem perder muitos graus de liberdade. No entanto, vale a pena considerar se o uso de uma variável fator envolveria um ajuste excessivo dos dados. Dado que existe uma relação bastante suave em seu gráfico, uma função linear simples ou quadrática obteria bons resultados sem ajuste excessivo.

Como você faz Participantuma inclinação aleatória, mas não intercepta a variável?

Você já colocou sua variável de resposta em uma escala de log usando uma função de link logarítmico, portanto, um efeito de interceptação para Participantestá dando um efeito multiplicativo na resposta. Se você der a isso uma inclinação aleatória interagindo com NumberContactsela, isso terá um efeito baseado em energia na resposta. Se você quiser isso, poderá obtê-lo com o (~ -1 + NumberContacts|Participant)que removerá a interceptação, mas adicionará uma inclinação com base no número de contatos.

Devo usar o Box-Cox para transformar meus dados? (por exemplo, lambda = 0,779)

$\lambda$

Devo incluir pesos para variação?

Comece examinando seu gráfico residual para ver se há evidências de heterocedasticidade. Com base nas plotagens que você já incluiu, parece-me que isso não é um problema; portanto, você não precisa adicionar pesos para a variação. Em caso de dúvida, você pode adicionar pesos usando uma função linear simples e, em seguida, executar um teste estatístico para verificar se a inclinação da ponderação é plana. Isso equivaleria a um teste formal de heterocedasticidade, que daria a você algum backup para sua escolha.

Devo incluir a autocorrelação NumberContacts?

Se você já incluiu um termo de efeito aleatório para o participante, provavelmente seria uma má idéia adicionar um termo de correlação automática ao número de contatos. Sua experiência usa um dedo diferente para diferentes números de contatos, para que você não espere a autocorrelação para o caso em que já contabilizou o participante. Adicionar um termo de autocorrelação além do efeito do participante significaria que você acha que existe uma dependência condicional entre o resultado de diferentes dedos, com base no número de contatos, mesmo para um determinado participante.

$^\dagger$

De fato, é razoável argumentar que as medidas tiradas de um participante não são independentes daquelas tiradas de outro participante. Por exemplo, algumas pessoas tendem a pressionar o dedo com mais (ou menos) força, o que afetaria todas as suas medições em cada número de contatos.

Portanto, a ANOVA de medidas repetidas de duas vias seria um modelo aceitável para aplicar neste caso.

Como alternativa, também se poderia aplicar um modelo de efeitos mistos participantcomo fator aleatório. Essa é uma solução mais avançada e sofisticada.