Bacteria picked up on fingers after multiple surface contacts: non-normal data, repeated measures, crossed participants

Question

Intro

I have participants who are repeatedly touching contaminated surfaces with E. coli in two conditions (A=wearing gloves, B=no gloves). I want to know if there's a difference between the amount of bacteria on their fingertips with and without gloves but also between the number of contacts. Both factors are within-participant.

Experimental Method:

Participants (n=35) touch each square once with the same finger for a maximum of 8 contacts (see figure a).

I then swab the finger of the participant and measure the bacteria on the fingertip after each contact. They then use a new finger to touch a different number of surfaces and so on from 1 to 8 contacts (see figure b).

Here is the real data: real data

The data is non-normal so see marginal distribution of bacteria|NumberContacts below. x=bacteria. Each facet is a different number of contacts.

MODEL

Trying from lme4::glmer based on amoeba's suggestions using Gamma(link="log") and polynomial for NumberContacts:

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

NB. Gamma(link="inverse") won't run saying PIRLS step-halving failed to reduce deviance.

Results:

Fitted vs residuals for cfug

qqp(resid(cfug))

Question:

Is my glmer model properly defined to incorporate the random effects of each participant and the fact that everyone does both experiment A followed by experiment B?

Addition:

Autocorrelation seems to exist between participants. This is probably because they weren't tested on the same day and the flask of bacteria grows and declines over time. Does it matter?

acf(CFU,lag=35) shows an significant correlation between one participant and the next.

Answer 1

Some plots to explore the data

Below are eight, one for each number of surface contacts, x-y plots showing gloves versus no gloves.

Each individuals is plotted with a dot. The mean and variance and covariance are indicated with a red dot and the ellipse (Mahalanobis distance corresponding to 97.5% of the population).

You can see that the effects are only small in comparison to the spread of the population. The mean is higher for 'no gloves' and the mean changes a bit higher up for more surface contacts (which can be shown to be significant). But the effect is only little in size (overall a $\frac{1}{4}$ log reduction), and there are many individuals for who there is actually a higher bacteria count with the gloves.

The small correlation shows that there is indeed a random effect from the individuals (if there wasn't an effect from the person then there should be no correlation between the paired gloves and no gloves). But it is only a small effect and an individual may have different random effects for 'gloves' and 'no gloves' (e.g. for all different contact point the individual may have consistently higher/lower counts for 'gloves' than 'no gloves').

Below plot are separate plots for each of the 35 individuals. The idea of this plot is to see whether the behavior is homogeneous and also to see what kind of function seems suitable.

Note that the 'without gloves' is in red. In most of the cases the red line is higher, more bacteria for the cases 'without gloves'.

I believe that a linear plot should be sufficient to capture the trends here. The disadvantage of the quadratic plot is that the coefficients are gonna be more difficult to interpret (you are not gonna see directly whether the slope is positive or negative because both the linear term and the quadratic term have an influence on this).

But more importantly you see that the trends differ a lot among the different individuals and therefore it may be useful to add a random effect for not only the intercept, but also the slope of the individual.

Model

With the model below

• Each individual will get it's own curve fitted (random effects for linear coefficients).
• The model uses log-transformed data and fits with a regular (gaussian) linear model. In the comments amoeba mentioned that a log link is not relating to a lognormal distribution. But this is different. $y \sim N(\log(\mu),\sigma^2)$ is different from $\log(y) \sim N(\mu,\sigma^2)$
• Weights are applied because the data is heteroskedastic. The variation is more narrow towards the higher numbers. This is probably because the bacteria count has some ceiling and the variation is mostly due to failing transmission from the surface to the finger (= related to lower counts). See also in the 35 plots. There are mainly a few individuals for which the variation is much higher than the others. (we see also bigger tails, overdispersion, in the qq-plots)
• No intercept term is used and a 'contrast' term is added. This is done to make the coefficients easier to interpret.

.

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)

This gives

> 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

code to obtain plots

chemometrics::drawMahal function

# 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')
}

Answer 2

Firstly, good work on your graph; it gives a clear representation of the data, so you can already see the kind of pattern in the data based on the number of contacts and the use or absence of gloves.$^\dagger$ Looking at this graph, I would think you will get good results with a basic log-polynomial model, with random effects for the participants. The model you have chosen looks reasonable, but you might also consider adding a quadratic term for the number of contacts.

As to whether to use MASS:glmmPQL or lme4:glmer for your model, my understanding is that both these functions will fit the same model (so long as you set the model equation, distribution and link function the same) but they use different estimation methods to find the fit. I could be mistaken, but my understanding from the documentation is that glmmPQL uses penalised quasi-likelihood as described in Wolfinger and O'Connell (1993), whereas glmer uses Gauss-Hermite quadrature. If you are worried about it you can fit your model with both methods and check that they give the same coefficient estimates and that way you will have greater confidence that the fitting algorithm has converged to the true MLEs of the coefficients.

---

Should NumberContacts be a categorical factor?

This variable has a natural ordering that appears from your plots to have a smooth relationship with the response variable, so you could reasonably treat it as a numeric variable. If you were to include factor(NumberContacts) then you will not constrain its form and you will not lose many degrees-of-freedom. You could even use the interaction Gloves*factor(NumberContacts) without losing too many degrees-of-freedom. However, it is worth considering whether using a factor variable would involve over-fitting the data. Given that there is a fairly smooth relationship in your plot a simple linear function or quadratic would get good results without over-fitting.

---

How do you make Participant a random slope but not intercept variable?

You have already put your response variable on a log-scale by using a logarithmic link function, so an intercept effect for Participant is giving a multiplicative effect on the response. If you were to give this a random slope interacting with NumberContacts then it would have a power-based effect on the response. If you want this then you can get it with (~ -1 + NumberContacts|Participant) which will remove the intercept but add a slope based on the number of contacts.

---

Should I use Box-Cox for transforming my data? (e.g. lambda=0.779)

If in doubt, try fitting a model with this transformation and see how it compares to other models using appropriate goodness-of-fit statistics. If you are going to use this transformation it would be better to leave the parameter $\lambda$ as a free parameter and let it be estimated as part of your model, rather than pre-specifying a value.

---

Should I include weights for variance?

Start by looking at your residual plot to see if there is evidence of heteroscedasticity. Based on the plots you have already included it looks to me like this is not an issue, so you don't need to add any weights for the variance. If in doubt you could add weights using a simple linear function and then perform a statistical test to see if the slope of the weighting is flat. That would amount to a formal test of heteroscedasticity, which would give you some backup for your choice.

---

Should I include autocorrelation in NumberContacts?

If you have already included a random effect term for the participant then it would probably be a bad idea to add an auto-correlation term on the number of contacts. Your experiment uses a different finger for different numbers of contacts so you would not expect autocorrelation for the case where you have already accounted for the participant. Adding an autocorrelation term in addition to the participant effect would mean that you think there is a conditional dependence between the outcome of different fingers, based on the number of contacts, even for a given participant.

---

$^\dagger$ Your graph is good in showing the relationships, but you could improve it aesthetically by adding a title and subtitle information and giving it better axis labels. You could also simplify your legend by removing its title and changing 'Yes' to 'Gloves' and 'No' to 'No Gloves'.

Answer 3

Indeed, it is reasonable to argue that measurements taken from one participant are not independent from those taken from another participant. For example, some people might tend to press their finger with more (or less) force, which would affect all their measurements accross each number of contacts.

So the 2-way repeated-measures ANOVA would be an acceptable model to apply in this case.

Alternatively, one could also apply a mixed-effects model with participant as a random factor. This is a more advanced and a more sophisticated solution.