Bakteri mengambil dengan jari setelah beberapa kontak permukaan: data yang tidak normal, tindakan berulang, peserta silang

Intro

Saya memiliki peserta yang berulang kali menyentuh permukaan yang terkontaminasi dengan E. coli dalam dua kondisi ( A = memakai sarung tangan, B = tanpa sarung tangan). Saya ingin tahu apakah ada perbedaan antara jumlah bakteri di ujung jari mereka dengan dan tanpa sarung tangan tetapi juga antara jumlah kontak. Kedua faktor tersebut adalah di dalam partisipan.

Metode eksperimen:

Peserta (n = 35) menyentuh setiap kotak sekali dengan jari yang sama untuk maksimum 8 kontak (lihat gambar a).

Saya kemudian menyeka jari peserta dan mengukur bakteri di ujung jari setelah setiap kontak. Mereka kemudian menggunakan jari baru untuk menyentuh jumlah permukaan yang berbeda dan seterusnya dari 1 hingga 8 kontak (lihat gambar b).

Berikut adalah data nyata : data nyata

Data ini tidak normal jadi lihat distribusi marginal bakteri | NumberContacts di bawah ini. x = bakteri. Setiap sisi adalah jumlah kontak yang berbeda.

MODEL

Mencoba dari lme4 :: glmer berdasarkan saran amoeba menggunakan Gamma (tautan = "log") dan polinomial untuk NumberContacts:

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

NB. Gamma (tautan = "terbalik") tidak akan berjalan dengan mengatakan PIRLS separuh langkah gagal mengurangi penyimpangan.

Hasil:

Dipasang vs residu untuk cfug

qqp (resid (cfug))

Pertanyaan:

Apakah model glmer saya didefinisikan dengan benar untuk menggabungkan efek acak dari setiap peserta dan fakta bahwa setiap orang melakukan keduanya percobaan A diikuti oleh percobaan B ?

Tambahan:

Autokorelasi tampaknya ada di antara peserta. Ini mungkin karena mereka tidak diuji pada hari yang sama dan labu bakteri tumbuh dan menurun dari waktu ke waktu. Apakah itu penting?

ACF (CFU, lag = 35) menunjukkan korelasi yang signifikan antara satu peserta dan yang berikutnya.

Beberapa plot untuk mengeksplorasi data

Di bawah ini delapan, satu untuk setiap jumlah kontak permukaan, petak xy yang menunjukkan sarung tangan versus tanpa sarung tangan.

Setiap individu diplot dengan sebuah titik. Rerata dan varians dan kovarians ditunjukkan dengan titik merah dan elips (jarak Mahalanobis sesuai dengan 97,5% populasi).

$\frac{1}{4}$

Korelasi kecil menunjukkan bahwa memang ada efek acak dari individu (jika tidak ada efek dari orang tersebut maka tidak boleh ada korelasi antara sarung tangan berpasangan dan tidak ada sarung tangan). Tetapi ini hanya efek kecil dan seorang individu mungkin memiliki efek acak berbeda untuk 'sarung tangan' dan 'tidak ada sarung tangan' (misalnya untuk semua titik kontak yang berbeda, individu mungkin memiliki jumlah yang lebih tinggi / lebih rendah secara konsisten untuk 'sarung tangan' daripada 'tanpa sarung tangan') .

Plot di bawah ini adalah plot terpisah untuk masing-masing 35 individu. Gagasan plot ini adalah untuk melihat apakah perilaku itu homogen dan juga untuk melihat fungsi apa yang cocok.

Perhatikan bahwa 'tanpa sarung tangan' berwarna merah. Dalam sebagian besar kasus, garis merah lebih tinggi, lebih banyak bakteri untuk kasus 'tanpa sarung tangan'.

Saya percaya bahwa plot linier harus cukup untuk menangkap tren di sini. Kerugian dari plot kuadratik adalah bahwa koefisien akan lebih sulit untuk ditafsirkan (Anda tidak akan melihat secara langsung apakah kemiringan positif atau negatif karena istilah linear dan kuadratik memiliki pengaruh pada hal ini).

Tetapi yang lebih penting Anda melihat bahwa tren sangat berbeda di antara individu yang berbeda dan karena itu mungkin berguna untuk menambahkan efek acak untuk tidak hanya intersep, tetapi juga kemiringan individu.

Model

Dengan model di bawah ini

• Setiap individu akan mendapatkan kurva itu sendiri (efek acak untuk koefisien linier).
• $y \sim N(\log(\mu),\sigma^2)$$\log(y) \sim N(\mu,\sigma^2)$
• Bobot diterapkan karena datanya heteroskedastik. Variasi lebih sempit ke arah angka yang lebih tinggi. Ini mungkin karena jumlah bakteri memiliki langit-langit dan variasi ini sebagian besar disebabkan oleh kegagalan transmisi dari permukaan ke jari (= terkait dengan jumlah yang lebih rendah). Lihat juga di 35 plot. Ada beberapa individu yang variasinya jauh lebih tinggi daripada yang lain. (kami juga melihat ekor yang lebih besar, penyebaran berlebihan, dalam plot qq)
• Tidak ada istilah intersepsi yang digunakan dan istilah 'kontras' ditambahkan. Ini dilakukan untuk membuat koefisien lebih mudah untuk ditafsirkan.

.

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)

Ini memberi

> 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

kode untuk mendapatkan plot

chemometrics :: function 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 plot

#### 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 x 4 plot

#### 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$

Mengenai apakah akan menggunakan MASS:glmmPQLatau lme4:glmeruntuk model Anda, pemahaman saya adalah bahwa kedua fungsi ini akan cocok dengan model yang sama (selama Anda mengatur persamaan model, distribusi dan fungsi tautan sama), tetapi mereka menggunakan metode estimasi yang berbeda untuk menemukan kecocokan. Saya bisa saja salah, tetapi pemahaman saya dari dokumentasi adalah yang glmmPQLmenggunakan kuasi-kemungkinan hukuman seperti yang dijelaskan dalam Wolfinger dan O'Connell (1993) , sedangkan glmermenggunakan quadrature Gauss-Hermite. Jika Anda khawatir tentang hal itu, Anda dapat mencocokkan model Anda dengan kedua metode dan memeriksa apakah mereka memberikan estimasi koefisien yang sama dan dengan cara itu Anda akan memiliki keyakinan yang lebih besar bahwa algoritma pemasangan telah konvergen ke MLEs sebenarnya dari koefisien.

Haruskah NumberContactsmenjadi faktor kategorikal?

Variabel ini memiliki urutan alami yang muncul dari plot Anda untuk memiliki hubungan yang lancar dengan variabel respons, sehingga Anda dapat memperlakukannya sebagai variabel numerik. Jika Anda memasukkannya factor(NumberContacts)maka Anda tidak akan membatasi bentuknya dan Anda tidak akan kehilangan banyak derajat kebebasan. Anda bahkan bisa menggunakan interaksi Gloves*factor(NumberContacts)tanpa kehilangan terlalu banyak derajat kebebasan. Namun, perlu dipertimbangkan apakah menggunakan variabel faktor akan melibatkan pemasangan data yang berlebihan. Mengingat bahwa ada hubungan yang cukup lancar dalam plot Anda, fungsi linear sederhana atau kuadratik akan mendapatkan hasil yang baik tanpa pemasangan berlebihan.

Bagaimana Anda membuat Participantkemiringan acak tetapi tidak mencegat variabel?

Anda telah meletakkan variabel respons Anda pada skala log dengan menggunakan fungsi tautan logaritmik, sehingga efek intersep untuk Participantmemberi efek multiplikatif pada respons. Jika Anda memberi ini kemiringan acak yang berinteraksi dengannya NumberContactsmaka itu akan memiliki efek berbasis kekuatan pada respons. Jika Anda menginginkan ini maka Anda bisa mendapatkannya dengan (~ -1 + NumberContacts|Participant)yang akan menghapus intersepsi tetapi menambahkan kemiringan berdasarkan jumlah kontak.

Haruskah saya menggunakan Box-Cox untuk mengubah data saya? (misalnya lambda = 0,779)

$\lambda$

Haruskah saya memasukkan bobot untuk varian?

Mulailah dengan melihat plot residu Anda untuk melihat apakah ada bukti heteroskedastisitas. Berdasarkan plot yang sudah Anda sertakan, tampaknya bagi saya ini bukan masalah, jadi Anda tidak perlu menambahkan bobot untuk varians. Jika ragu Anda bisa menambahkan bobot menggunakan fungsi linier sederhana dan kemudian melakukan uji statistik untuk melihat apakah kemiringan bobotnya rata. Itu akan menjadi tes formal heteroskedastisitas, yang akan memberi Anda cadangan untuk pilihan Anda.

Haruskah saya menyertakan autokorelasi NumberContacts?

Jika Anda sudah memasukkan istilah efek acak untuk peserta, maka mungkin ide buruk untuk menambahkan istilah korelasi otomatis pada jumlah kontak. Percobaan Anda menggunakan jari yang berbeda untuk jumlah kontak yang berbeda sehingga Anda tidak akan mengharapkan autokorelasi untuk kasus di mana Anda telah memperhitungkan peserta. Menambahkan istilah autokorelasi di samping efek peserta akan berarti bahwa Anda berpikir ada ketergantungan bersyarat antara hasil jari yang berbeda, berdasarkan jumlah kontak, bahkan untuk peserta yang diberikan.

$^\dagger$

Memang, masuk akal untuk berpendapat bahwa pengukuran yang diambil dari satu peserta tidak independen dari yang diambil dari peserta lain. Sebagai contoh, beberapa orang mungkin cenderung menekan jari mereka dengan kekuatan lebih (atau kurang), yang akan memengaruhi semua pengukuran mereka di setiap jumlah kontak.

Jadi ANOVA 2 arah yang diulang akan menjadi model yang dapat diterima untuk diterapkan dalam kasus ini.

Atau, kita juga bisa menerapkan model efek campuran dengan participantsebagai faktor acak. Ini adalah solusi yang lebih maju dan lebih canggih.