Niektóre wykresy do eksploracji danych
Poniżej znajduje się osiem, po jednym dla każdej liczby styków powierzchniowych, wykresy xy pokazujące rękawice w porównaniu z brakiem rękawic.
Każda osoba jest wykreślona kropką. Średnia, wariancja i kowariancja są oznaczone czerwoną kropką i elipsą (odległość Mahalanobisa odpowiadająca 97,5% populacji).
Widać, że efekty są niewielkie w porównaniu do rozprzestrzeniania się populacji. Średnia jest wyższa dla „bez rękawiczek”, a średnia zmienia się nieco wyżej w przypadku większej liczby kontaktów na powierzchni (co można wykazać jako znaczące). Ale efekt jest tylko trochę w rozmiarze (ogólnie obniżenie log), i istnieje wiele osób, bo kto tam jest w rzeczywistości wyższe bakterie liczyć z rękawic.14
Mała korelacja pokazuje, że rzeczywiście istnieje efekt losowy od poszczególnych osób (jeśli nie wystąpiłby efekt od tej osoby, wówczas nie powinno być żadnej korelacji między sparowanymi rękawiczkami i bez rękawiczek). Jest to jednak tylko niewielki efekt, a osoba może mieć różne losowe skutki dla „rękawiczek” i „bez rękawiczek” (np. Dla wszystkich różnych punktów kontaktowych dana osoba może mieć stale wyższą / niższą liczbę dla „rękawiczek” niż „bez rękawiczek”) .
Poniżej wykresu są osobne wykresy dla każdej z 35 osobników. Ideą tego wykresu jest sprawdzenie, czy zachowanie jest jednorodne, a także sprawdzenie, jaki rodzaj funkcji wydaje się odpowiedni.
Pamiętaj, że „bez rękawiczek” ma kolor czerwony. W większości przypadków czerwona linia jest wyższa, więcej bakterii w przypadkach „bez rękawiczek”.
Uważam, że liniowa fabuła powinna wystarczyć do uchwycenia trendów tutaj. Wadą wykresu kwadratowego jest to, że współczynniki będą trudniejsze do interpretacji (nie zobaczysz bezpośrednio, czy nachylenie jest dodatnie czy ujemne, ponieważ wpływ na to mają zarówno element liniowy, jak i kwadratowy).
Ale co ważniejsze, widać, że trendy różnią się znacznie między poszczególnymi osobami, dlatego może być użyteczny dodanie losowego efektu nie tylko dla przechwytywania, ale także nachylenia jednostki.
Model
Z poniższym modelem
- Każda osoba otrzyma własną krzywą (efekty losowe dla współczynników liniowych).
- Model wykorzystuje dane przekształcone w log i pasuje do zwykłego (gaussowskiego) modelu liniowego. W komentarzach amoeba wspomniała, że link do dziennika nie odnosi się do rozkładu logarytmicznego. Ale to jest inne.y∼ N.( log( Μ ) ,σ2)) jest różny od log( y) ∼ N( μ ,σ2))
- Wagi są stosowane, ponieważ dane są heteroskedastyczne. Różnica jest węższa w stosunku do wyższych liczb. Wynika to prawdopodobnie z tego, że liczba bakterii ma pewien pułap, a zmienność wynika głównie z nieudanej transmisji z powierzchni na palec (= związane z mniejszą liczbą). Zobacz także na 35 działkach. Istnieje głównie kilka osobników, dla których zróżnicowanie jest znacznie wyższe niż u innych. (widzimy także większe ogony, nadmierną dyspersję, na wykresach qq)
- Nie stosuje się terminu przechwytującego i dodaje się termin „kontrastowy”. Ma to na celu ułatwienie interpretacji współczynników.
.
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)
To daje
> 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
kod do uzyskania wykresów
chemometrics :: funkcja 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()
}
Działka 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))
}
Działka 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')
}
NumberContactsjako czynnika liczbowego i dołączyć kwadratowy / sześcienny termin wielomianowy. Lub spójrz na Uogólnione mieszane modele dodatków.