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

9

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). a) contatos dos dedos com 8 superfícies, b) UFC nos dedos após cada contato da superfície

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.

insira a descrição da imagem aqui

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 insira a descrição da imagem aqui

qqp (resid (cfug))

insira a descrição da imagem aqui

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.

insira a descrição da imagem aqui

HCAI
fonte
11
Você pode usar NumberContactscomo um fator numérico e incluir termos polinomiais quadráticos / cúbicos. Ou consulte Modelos mistos aditivos generalizados.
Ameba
11
@amoeba Obrigado por sua ajuda. Todos os participantes fizeram B (sem luva) seguido de A (com luva). Você acha que existem outros problemas fundamentais com a análise? Nesse caso, estou aberto a outras respostas.
HCAI
11
Nesse caso, você pode incluir o efeito aleatório da luva. Além disso, não entendo por que você remove a interceptação aleatória e por que não inclui todo o polinômio de 2º grau na parte aleatória. E você pode ter uma interação luva * num. Então, por que não CFU ~ Gloves * poly(NumberContacts,2) + (Gloves * poly(NumberContacts,2) | Participant)ou algo assim.
Ameba
11
Ah, eu entendo sobre a interceptação, mas você também precisaria suprimir a interceptação fixa. Além disso, para zero contatos, você deve ter zero CFU, mas com o link de log isso não faz sentido. E você não tem nem perto de zero CFU em um contato. Portanto, eu não suprimia a interceptação. Não convergindo não é bom, tentar remover a interação da parte aleatória: CFU ~ Gloves * poly(NumberContacts,2) + (Gloves + poly(NumberContacts,2) | Participant), ou talvez remover as luvas de lá CFU ~ Gloves * poly(NumberContacts,2) + (poly(NumberContacts,2) | Participant)...
ameba
11
Eu acho que Gloves * poly(NumberContacts,2) + (poly(NumberContacts,2) | Participant)é um modelo bastante decente.
Ameba

Respostas:

6

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.14

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') .

xy parcelas de com e 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.

parcelas para cada 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 )yN(log(μ),σ2)log(y)N(μ,σ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

qqplot

resíduos

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')
}
Sextus Empiricus
fonte
Muito obrigado, Martijn, você explicou as coisas com tanta clareza. Surpreendente! Como a recompensa terminou antes que eu pudesse atribuí-la, gostaria muito de lhe oferecer uma quantia separada (veremos como fazer isso agora). Tenho algumas perguntas: Em primeiro lugar, transformar os dados parece ter escolas de pensamento: alguns concordam e outros discordam veementemente. Por que está tudo bem aqui? Em segundo lugar, por que remover a interceptação aleatória facilita a interpretação dos coeficientes?
HCAI
(2) Eu acho que a transformação está correta quando você poderia argumentar que existe um processo que torna a transformação lógica (na verdade, a transformação relutante porque faz com que os resultados pareçam bons pode ser vista como manipulação de dados e deturpação de resultados, além de não obter o subjacente )
Sextus Empiricus
Vejo @Martijn, pelo menos em biologia, a transformação por log10 é comum para bactérias. Fico feliz em dar a recompensa, você merece. Você se importaria de elaborar um pouco o porquê de usar esse "termo de contraste", por favor?
HCAI
11
Com relação ao contraste Consulte aqui stats.stackexchange.com/a/308644/164061 Você tem a liberdade de mudar o termo de interceptação. Uma maneira possivelmente útil é definir a interceptação entre as duas categorias e deixar que o efeito seja a diferença entre os dois efeitos (um será negativo e outro positivo) em relação ao termo médio de interceptação. (não que eu tinha para adicionar uma variável para isso)
Sexto Empírico
11
Idealmente, você teria os tratamentos distribuídos aleatoriamente ao longo do tempo, para que quaisquer efeitos possíveis devido a variações no tempo se nivelassem. Mas, na verdade, não vejo tanta autocorrelação. Você quer dizer saltos como no participante 5, entre 5 e 6, número de contatos após o qual a linha está estável novamente? Eu acho que essas não são tão ruins e, no máximo, adicionam ruído, mas não interferem no seu método (exceto fazendo com que o sinal / ruído seja baixo). Você pode ter mais certeza quando não vê uma mudança sistemática ao longo do tempo. Se você processou os participantes em ordem, pode plotar a média de CFU ao longo do tempo.
Sextus Empiricus
2

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)

λ


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.


Ben - Restabelecer Monica
fonte
Obrigado, essa é uma resposta incrível! No final, tentei Gamma (link = "log") e o glmer converge sem reclamar, viva! glmer (CFU ~ Luvas + poli (NumberContacts, 2) + (-1 + NumberContacts | Participante), dados = na.omit (subconjunto (K, CFU <4.5e5 & Surface == "P")), família = Gamma ( link = "log")). QQplot acho que está OK (nada fora dos ICs), mas os redimensionados estão ajustados (veja a foto adicionada após o comentário ser publicado, caso não corresponda). Devo me preocupar muito com isso?
HCAI
11
O enredo de QQ parece bom para mim. Além disso, lembre-se de que em um GLM os resíduos de Pearson não seguem necessariamente uma distribuição normal. Parece que você tem uma boa análise.
Ben - Restabelece Monica
1

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.

Mihael
fonte
Obrigado Mihael, você está absolutamente certo sobre a pressão. Hmm, eu estava lendo sobre o modelo de efeitos mistos aqui rcompanion.org/handbook/I_09.html, mas não tenho certeza sobre as interações e os fatores aninhados. Meus fatores estão aninhados?
HCAI
Devo também salientar que os dados normalmente não são distribuídos para cada contato, portanto, analisamos a modelagem de Quase-Probabilidade Penalizada (PQL): ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/… . Você acha que essa é uma boa escolha?
HCAI