Бактерии, обнаруженные на пальцах после нескольких поверхностных контактов: ненормальные данные, повторные измерения, скрещенные участники

вступление

У меня есть участники, которые неоднократно касаются загрязненных поверхностей кишечной палочкой в ​​двух состояниях ( A = ношение перчаток, B = отсутствие перчаток). Я хочу знать, есть ли разница между количеством бактерий на кончиках пальцев в перчатках и без них, а также между количеством контактов. Оба фактора находятся внутри участника.

Экспериментальный метод:

Участники (n = 35) касаются каждого квадрата один раз одним и тем же пальцем максимум для 8 контактов (см. Рисунок а).

Затем я протираю палец участника и измеряю бактерии на кончике пальца после каждого контакта. Затем они с помощью нового пальца касаются разного количества поверхностей и т. Д. От 1 до 8 контактов (см. Рисунок б).

Вот реальные данные: реальные данные

Данные не являются нормальными, поэтому смотрите предельное распределение бактерий | NumberContacts ниже. х = бактерия. Каждый аспект - это разное количество контактов.

МОДЕЛЬ

Попытка из lme4 :: glmer на основе предложений amoeba с использованием Gamma (link = "log") и полинома для NumberContacts:

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

NB. Гамма (link = "inverse") не будет работать, говоря, что разделение на два шага PIRLS не смогло уменьшить отклонение.

Результаты:

Приспособлено против остатков для cfug

qqp (кубовые остатки (cfug))

Вопрос:

Правильно ли определена моя модель блеска, чтобы учесть случайные эффекты каждого участника и тот факт, что каждый проводит эксперимент A с последующим экспериментом B ?

Дополнение:

Автокорреляция, кажется, существует между участниками. Вероятно, это связано с тем, что они не были проверены в один и тот же день, а колба бактерий растет и со временем уменьшается. Это имеет значение?

acf (CFU, лаг = 35) показывает значительную корреляцию между одним участником и следующим.

Некоторые участки для изучения данных

Ниже приведены восемь, по одному для каждого числа поверхностных контактов, на графиках xy показаны перчатки, а не перчатки.

Каждый человек нанесен на карту с точкой. Среднее значение, дисперсия и ковариация обозначены красной точкой и эллипсом (расстояние Махаланобиса соответствует 97,5% населения).

Вы можете видеть, что эффекты только незначительны по сравнению с распространением населения. Среднее значение выше для «без перчаток», и среднее значение немного выше для большего количества поверхностных контактов (что может быть показано как существенное). Но эффект только небольшой по размеру (в целом$\frac{1}{4}$ сокращение журнала), и есть много людей , для тех , кто есть на самом деле рассчитывать более высокие бактерии с перчатками.

Небольшая корреляция показывает, что у людей действительно есть случайный эффект (если у человека не было эффекта, то не должно быть корреляции между парными перчатками и перчатками). Но это лишь небольшой эффект, и у индивидуума могут быть разные случайные эффекты для «перчаток» и «без перчаток» (например, для всех разных точек контакта у индивидуума могут быть постоянно более высокие / более низкие показатели для «перчаток», чем «без перчаток») ,

Ниже приведены отдельные участки для каждого из 35 человек. Идея этого графика заключается в том, чтобы увидеть, является ли поведение однородным, а также посмотреть, какая функция кажется подходящей.

Обратите внимание, что «без перчаток» красного цвета. В большинстве случаев красная линия выше, больше бактерий для случаев «без перчаток».

Я считаю, что линейного сюжета должно быть достаточно, чтобы уловить здесь тенденции. Недостаток квадратичного графика состоит в том, что коэффициенты будет труднее интерпретировать (вы не увидите непосредственно, является ли наклон положительным или отрицательным, потому что на это влияют и линейный член, и квадратный член).

Но что еще более важно, вы видите, что тенденции сильно различаются между разными людьми, и поэтому может быть полезно добавить случайный эффект не только для перехвата, но и для склонности человека.

модель

С моделью ниже

• Каждый индивид получит свою собственную кривую (случайные эффекты для линейных коэффициентов).
• $y \sim N(\log(\mu),\sigma^2)$$\log(y) \sim N(\mu,\sigma^2)$
• Веса применяются потому, что данные гетероскедастичны. Вариация более узкая в сторону больших чисел. Вероятно, это связано с тем, что количество бактерий имеет некоторый потолок, и отклонения в основном связаны с отсутствием передачи с поверхности на палец (= связано с меньшим количеством). Смотрите также на 35 графиках. Есть в основном несколько человек, у которых разброс намного выше, чем у других. (мы видим также большие хвосты, избыточную дисперсию на qq-графиках)
• Термин «перехват» не используется, и добавляется термин «контраст». Это сделано для облегчения интерпретации коэффициентов.

,

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)

Это дает

> 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

код для получения участков

хемометрика :: функция 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 х 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))
}

2 х 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$ Глядя на этот график, я думаю, вы получите хорошие результаты с базовой лог-полиномиальной моделью со случайными эффектами для участников. Модель, которую вы выбрали, выглядит разумной, но вы также можете рассмотреть вопрос о добавлении квадратичного термина для числа контактов.

Что касается использования MASS:glmmPQLили lme4:glmerвашей модели, я понимаю, что обе эти функции будут соответствовать одной и той же модели (при условии, что вы задаете уравнение модели, распределение и функцию связи одинаково), но они используют разные методы оценки для нахождения соответствия. Я могу ошибаться, но, насколько я понимаю из документации, glmmPQLиспользуется штрафная квази-правдоподобие, как описано в Wolfinger and O'Connell (1993) , тогда как glmerиспользуется квадратура Гаусса-Эрмита. Если вы беспокоитесь об этом, вы можете согласовать свою модель с обоими методами и убедиться, что они дают одинаковые оценки коэффициентов, и таким образом у вас будет большая уверенность в том, что алгоритм подбора сходится к истинным MLE коэффициентов.

Должен NumberContactsбыть категоричный фактор?

Эта переменная имеет естественный порядок, который, как видно из ваших графиков, имеет гладкую связь с переменной отклика, поэтому вы можете разумно рассматривать ее как числовую переменную. Если вы factor(NumberContacts)включите его, вы не будете ограничивать его форму и не потеряете много степеней свободы. Вы даже можете использовать взаимодействие, Gloves*factor(NumberContacts)не теряя слишком много степеней свободы. Тем не менее, стоит подумать, не приведет ли использование факторной переменной к переопределению данных. Учитывая, что в вашем графике есть довольно гладкие отношения, простая линейная или квадратичная функция могла бы получить хорошие результаты без чрезмерной подгонки.

Как вы делаете Participantслучайный наклон, но не переменную перехвата?

Вы уже поместили свою переменную ответа в логарифмическую шкалу, используя функцию логарифмической связи, поэтому эффект перехвата для Participantдает мультипликативный эффект на ответ. Если бы вы дали этому случайному наклону взаимодействовать с NumberContactsним, это бы имело эффект силы на ответ. Если вы хотите это, то вы можете получить его, (~ -1 + NumberContacts|Participant)который удалит перехват, но добавит наклон в зависимости от количества контактов.

Должен ли я использовать Box-Cox для преобразования моих данных? (например, лямбда = 0,779)

$\lambda$

Должен ли я включить веса для дисперсии?

Начните с просмотра вашего остаточного участка, чтобы увидеть, есть ли признаки гетероскедастичности. Исходя из графиков, которые вы уже включили, мне кажется, что это не проблема, поэтому вам не нужно добавлять весовые коэффициенты для дисперсии. Если вы сомневаетесь, вы можете добавить веса, используя простую линейную функцию, а затем выполнить статистический тест, чтобы увидеть, является ли наклон веса плоским. Это будет формальным тестом гетероскедастичности, который даст вам некоторую резервную копию для вашего выбора.

Должен ли я включить автокорреляцию в NumberContacts?

Если вы уже включили термин случайного эффекта для участника, то было бы неплохо добавить термин автокорреляции для числа контактов. В вашем эксперименте используется другой палец для разного количества контактов, поэтому вы не ожидаете автокорреляции для случая, когда вы уже учли участника. Добавление термина автокорреляции в дополнение к эффекту участника будет означать, что вы думаете, что между результатами разных пальцев существует условная зависимость, основанная на количестве контактов, даже для данного участника.

$^\dagger$

Действительно, разумно утверждать, что измерения, взятые у одного участника, не являются независимыми от измерений, проведенных у другого участника. Например, некоторые люди могут нажимать на палец с большей (или меньшей) силой, что повлияет на все их измерения по каждому количеству контактов.

Таким образом, ANOVA с двусторонним повторным измерением будет приемлемой моделью для применения в этом случае.

В качестве альтернативы можно также применить модель смешанных эффектов со participantслучайным фактором. Это более продвинутое и более сложное решение.