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

9

вступление

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

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

Участники (n = 35) касаются каждого квадрата один раз одним и тем же пальцем максимум для 8 контактов (см. Рисунок а). а) пальцевые контакты с 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) показывает значительную корреляцию между одним участником и следующим.

введите описание изображения здесь

HCAI
источник
1
Вы можете использовать NumberContactsв качестве числового множителя и включить квадратичные / кубические полиномиальные члены. Или посмотрите на Обобщенные Аддитивные Смешанные Модели.
амеба
1
@amoeba Спасибо за вашу помощь. Все участники сделали B (без перчаток), а затем A (в перчатках). Как вы думаете, есть другие фундаментальные проблемы с анализом? Если так, я открыт для любых дальнейших ответов.
HCAI
1
Если это так, то вы можете включить случайный эффект перчатки. Кроме того, я не понимаю, почему вы удаляете случайный перехват и почему вы не включаете весь полином 2-й степени в случайную часть. И вы можете иметь перчатки * Num взаимодействия. Так почему нет CFU ~ Gloves * poly(NumberContacts,2) + (Gloves * poly(NumberContacts,2) | Participant)или что-то в этом роде.
амеба
1
О, я понимаю о перехвате, но тогда вам также придется подавить фиксированный перехват. Кроме того, для нулевых контактов у вас должен быть нулевой КОЕ, но с лог-линком это не имеет смысла. И у вас нигде рядом нет нуля КОЕ на 1 контакт. Так что я бы не подавлял перехват. Не сходиться не хорошо, попробуйте удалить взаимодействие со случайной частью: CFU ~ Gloves * poly(NumberContacts,2) + (Gloves + poly(NumberContacts,2) | Participant)или, может быть, удалить Перчатки оттуда CFU ~ Gloves * poly(NumberContacts,2) + (poly(NumberContacts,2) | Participant)...
amoeba
1
Я думаю Gloves * poly(NumberContacts,2) + (poly(NumberContacts,2) | Participant), это довольно приличная модель.
амеба

Ответы:

6

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

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

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

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

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

xy участки с перчатками и без

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

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

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

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

участки для каждого человека

модель

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

  • Каждый индивид получит свою собственную кривую (случайные эффекты для линейных коэффициентов).
  • Y~N(журнал(μ),σ2)журнал(Y)~N(μ,σ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

qqplot

невязки

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

хемометрика :: функция 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')
}
Секст Эмпирик
источник
Большое вам спасибо, Мартин, вы так ясно все объяснили. Удивительно! Поскольку награда закончилась до того, как я смог назначить ее, я бы очень хотел предложить вам отдельную сумму (сейчас я посмотрю, как это сделать). У меня действительно есть несколько вопросов: во-первых, преобразование данных, кажется, вызывает сомнения: некоторые согласны, а некоторые категорически не согласны. Почему здесь хорошо? Во-вторых, почему удаление случайного пересечения облегчает интерпретацию коэффициентов?
HCAI
(2) Я полагаю, что преобразование в порядке, когда вы можете утверждать, что существует процесс, который делает преобразование логичным (действительно, преобразование неохотно, потому что оно делает результаты выглядящими хорошо, можно рассматривать как манипулирование данными и искажение результатов, а также отсутствие получения базового результата). модель)
Sextus Empiricus
Я вижу @Martijn, по крайней мере, в биологии трансформация с помощью log10 обычна для бактерий. Я рад отдать награду, вы это заслужили. Не могли бы вы немного рассказать о том, почему вы используете этот «контрастный термин», пожалуйста?
HCAI
1
Относительно контраста Смотрите здесь stats.stackexchange.com/a/308644/164061 У вас есть свобода перемещения термина перехвата. Один, возможно, полезный способ - установить перехват между двумя категориями и позволить эффекту быть разницей между двумя эффектами (один будет отрицательным, а другой положительным) относительно этого среднего члена перехвата. (не то чтобы я должен был добавить переменную для этого)
Sextus Empiricus
1
В идеале, вы должны иметь процедуры, случайно распределенные по времени, чтобы любые возможные эффекты из-за изменений во времени нивелировались. Но на самом деле я не вижу такой большой автокорреляции. Вы имеете в виду такие скачки, как у участника 5 между 5 и 6 числом контактов, после которого линия снова стабильна? Я думаю, что это не так уж плохо и в большинстве случаев добавляет шум, но не мешает вашему методу (за исключением того, что сигнал / шум низкий). Вы можете быть более уверены, когда не видите систематических изменений во времени. Если вы обработали участников по порядку, вы могли бы построить их среднее КОЕ со временем.
Секст Эмпирик
2

Глядя на этот график, я думаю, вы получите хорошие результаты с базовой лог-полиномиальной моделью со случайными эффектами для участников. Модель, которую вы выбрали, выглядит разумной, но вы также можете рассмотреть вопрос о добавлении квадратичного термина для числа контактов.

Что касается использования 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)

λ


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

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


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

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


Бен - Восстановить Монику
источник
Спасибо, это потрясающий ответ! В итоге я попробовал Gamma (link = "log") и блеск сходится без нареканий, ура! glmer (CFU ~ Перчатки + поли (NumberContacts, 2) + (-1 + NumberContacts | Участник), data = na.omit (подмножество (K, CFU <4.5e5 & Surface == "P")), семейство = Гамма ( ссылка = "журнал")). QQplot, я думаю, что все в порядке (ничего кроме CI), но приспособленные против остаточных лиц воронки (см. Добавленную картинку, добавленную после того, как этот комментарий опубликован, если он не совпадает). Должен ли я беспокоиться об этом слишком много?
HCAI
1
QQ сюжет выглядит хорошо для меня. Кроме того, помните, что в GLM остатки Пирсона не обязательно следуют нормальному распределению. Похоже, у вас есть хороший анализ.
Бен - Восстановить Монику
1

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

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

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

Mihael
источник
Спасибо, Майкл, ты абсолютно прав насчет давления. Хм, я читал здесь о модели смешанных эффектов rcompanion.org/handbook/I_09.html, но не уверен насчет взаимодействий и вложенных факторов. Являются ли мои факторы вложенными?
HCAI
Я также должен отметить, что данные обычно не распределяются по каждому контакту, поэтому обратились к моделированию штрафованных квази-правдоподобий (PQL): ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/… . Как вы думаете, это хороший выбор?
HCAI