R
多次表面接觸後手指上的細菌：非正常數據、重複測量、交叉參與者
介紹
我有參與者在兩種情況下反復接觸大腸桿菌污染的表面（A = 戴手套，B = 不戴手套）。我想知道戴手套和不戴手套時指尖上的細菌數量是否存在差異，以及接觸次數之間是否存在差異。這兩個因素都在參與者內部。
實驗方法：
參與者（n=35）用同一根手指觸摸每個方格一次，最多接觸 8 次（見圖 a）。
然後我用棉籤擦拭參與者的手指，並在每次接觸後測量指尖上的細菌。然後他們使用新手指觸摸不同數量的表面，依此類推，從 1 到 8 個觸點（見圖 b）。
這是真實數據：真實數據
數據是非正態的，因此請參閱下面的細菌邊緣分佈|NumberContacts。x=細菌。每個方面是不同數量的聯繫人。
模型
根據變形蟲的建議從lme4::glmer嘗試使用 Gamma(link=“log”) 和 NumberContacts 的多項式：
cfug<-glmer(CFU ~ Gloves + poly(NumberContacts,2) + (-1+NumberContacts|Participant), data=(K,CFU<4E5), family=Gamma(link="log") ) plot(cfug)注意。Gamma(link=“inverse”) 不會運行說 PIRLS 減半未能減少偏差。
結果：
cfug的擬合與殘差
qqp(resid(cfug))
題：
我的 glmer 模型是否正確定義以包含每個參與者的隨機效應以及每個人都做實驗A和實驗B的事實？
添加：
參與者之間似乎存在自相關。這可能是因為它們沒有在同一天進行測試，並且細菌瓶會隨著時間的推移而生長和減少。有關係嗎？
acf(CFU,lag=35) 顯示了一個參與者和下一個參與者之間的顯著相關性。
一些用於探索數據的圖
下面是八個，每個表面接觸數量一個，xy 圖顯示手套與不戴手套。
每個人都用一個點繪製。均值、方差和協方差用紅點和橢圓表示（馬氏距離對應於 97.5% 的總體）。
您可以看到，與人口分佈相比，影響很小。“不戴手套”的平均值更高，而更多表面接觸的平均值會更高（這可以證明是顯著的）。但效果只是很小（總體而言log 減少），並且有很多人實際上戴手套的細菌數量更高。
小的相關性表明確實存在來自個人的隨機效應（如果沒有來自人的影響，那麼配對手套和不戴手套之間應該沒有相關性）。但這只是一個很小的影響，個人可能對“手套”和“不戴手套”有不同的隨機影響（例如，對於所有不同的接觸點，個人對“手套”的計數可能始終高於/低於“不戴手套”） .
下圖是 35 個人中每個人的單獨圖。這個圖的想法是看行為是否是同質的，也看什麼樣的函數看起來合適。
請注意，“不戴手套”是紅色的。在大多數情況下，紅線越高，“不戴手套”的情況下細菌越多。
我相信線性圖應該足以捕捉這裡的趨勢。二次圖的缺點是係數將更難以解釋（您不會直接看到斜率是正數還是負數，因為線性項和二次項都會對此產生影響）。
但更重要的是，您會看到不同個體之間的趨勢差異很大，因此為截距和個體斜率添加隨機效應可能很有用。
模型
使用下面的模型
- 每個人都會得到自己的曲線擬合（線性係數的隨機效應）。
- 該模型使用對數轉換的數據並適合常規（高斯）線性模型。在評論中變形蟲提到日誌鏈接與對數正態分佈無關。但這是不同的。不同於
- 應用權重是因為數據是異方差的。數字越大，變化越窄。這可能是因為細菌計數有一定的上限，而變化主要是由於從表面到手指的傳播失敗（= 與較低的計數有關）。另見 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 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)) }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') }
