0
$\begingroup$

I have a small data set looking at the level of a particular protein (heat shock protein 70) in individuals exposed to different thermal treatments. There are 5 different thermal treatments with 10 individuals in each treatment. I'm interested to see if the thermal treatment affects the level of HSP70. Protein level is measured relative to a standard, so values vary from 0 to 1.2 (where 1 is the protein level of the standard).

In 2 of the 5 treatments, all individuals had no expression (i.e. 0 - it is likely that this is not a true zero, but rather below the detection of the equipment), in the other 3 treatments all individuals had varying levels of this protein. I had originally intended to analyze this as a one-way ANOVA, but with such a high amount of zeros (40%) in a small data set, it obviously doesn't pass normality or homogeneity of variances, even with various transformations. The non-parametric version of the one-way ANOVA (Krustal-Wallis) still requires homogeneity of variance, so that isn't appropriate either.

My dependent variable is (semi?) continuous so after googling for a while I came across zero-inflated semi continuous models, in particular a two-part or 'hurdle' model that might be appropriate. Is this overkill for such a small data set with one fixed factor?

How do I analyze this data? Should I use a different model? How do I do the zero-inflated semi continuous model in R (most zero-inflated models are associated with count data)?

$\endgroup$
1
  • $\begingroup$ Do you have a reasonable estimate of the detection limit? $\endgroup$
    – EdM
    Commented Apr 11, 2022 at 20:06

1 Answer 1

2
$\begingroup$

If you are measuring a continuous outcome with a lower limit of detection, it's best to think about the undetectable levels as left censored at the detection limit. That is, there are underlying values you can't detect because they are below the limit. This page is probably the most thorough discussion on this site.

The classic way to deal with a lower detection limit is tobit regression. There's a reasonable introduction with links on this page. You assume that there is some actual mean value with associated variance that is hidden for values below the detection limit. You use maximum likelihood to estimate the corresponding parameter values given that you can't observe anything below the limit. This page shows a way to implement this in R. Censored values are also an issue in survival analysis, and the methods used for tobit regressions are essentially survival models.

If you use a tobit model, be aware that its standard implementation assumes a Gaussian distribution of the values about the means, including the values that are hidden. That might lead to negative estimates of the undetected mean values, which wouldn't make much sense in your case. One way around that is to model the logarithms of the HSP70 values with a lognormal survival model. You model:

$$\log Y = X'\beta + \sigma W $$

where $Y$ is your outcome measure, $X$ is the vector of predictors (treatments in your case), $\beta$ is the corresponding vector of regression coefficients, $\sigma$ is a scale factor, and $W$ is a standard normal distribution. This is appropriate if, as is often the case in biochemical studies, the error is proportional to the actual value. The regression coefficients here represent fractional changes in the outcome.

Here's an example. Make up some data with 5 treatment groups and 10 observations each. Two groups have exactly 0 values for all observations and the other 3 groups have successively higher mean values of 0.4, 0.8, and 1.2 with a 25% coefficient of variation (SD/mean). We assume that the lower limit of detection is 0.2 units, so the observations with 0 values are given values of 0.2 but are noted as censored (observed = 0):

set.seed(555)
dfTobit <- data.frame(tx=rep(LETTERS[1:5],times=rep(10,5)),val=c(rep(0.2,20),rnorm(10,0.4,.1),rnorm(10,0.8,.2),rnorm(10,1.2,.3)),observed=c(rep(0,20),rep(1,30)))

Load the survival package and build a lognormal survival model. Specify that these are left-censored observations.

library(survival)
tobFit <- survreg(Surv(val,observed,type="left")~tx,data=dfTobit,dist="lognormal")

Calculate the residuals, which should match the assumed standard Gaussian distribution in the lognormal model. (Take the difference between the observed $\log Y$ values and the model's corresponding linear predictor estimates, then divide by the estimated scale factor.)

tobFit$scale
# [1] 0.212059
tobResids <- (log(dfTobit$val)-predict(tobFit,type="lp"))/.212

The Gaussian fit isn't bad except for one low outlier in group C (value of 0.22 versus a mean of 0.4) close to the cutoff of 0.2:

qqnorm(tobResids[dfTobit$observed==1],bty="n")
qqline(tobResids[dfTobit$observed==1])

qqnorm plot of tobit residuals

Predictions for the 5 groups and their standard errors:

predict(tobFit,newdata=data.frame(tx=c("A","B","C","D","E")),se.fit=TRUE)
# $fit
#          1          2          3          4          5 
# 0.05036946 0.05036946 0.42162094 0.69095441 1.23609443 
# 
# $se.fit
#           1           2           3           4           5 
# 81.75465693 81.75466962  0.02827346  0.04633468  0.08289121 

The standard errors are enormous for the two groups with 0 observed values, but the estimated mean values are non-negative. The mean estimates for the other 3 groups are the geometric rather than the raw means, as the model worked with log-transformed values.

You can do then pairwise comparisons with multiple-comparison correction, for example with the emmeans package, shown here for the ratios of the 3 treatments with non-zero values:

library(emmeans)
pairs(emmeans(tobFit,~tx),exclude=1:2,reverse=TRUE,type="response")
#  contrast ratio    SE df null t.ratio p.value
#  D / C     1.64 0.155 44    1   5.209  <.0001
#  E / C     2.93 0.278 44    1  11.342  <.0001
#  E / D     1.79 0.170 44    1   6.133  <.0001
# 
# P value adjustment: tukey method for comparing a family of 3 estimates 
# Tests are performed on the log scale 

The UCLA web page linked above shows how to do this type of analysis with the tobit() function of the VGAM package. I think that package also allows you to model the mean and the error term separately, for example if you wanted coefficients to provide additive effects while the errors were still proportional to the values. Apply your knowledge of the subject matter to decide the most appropriate distributions.

Other approaches

All that said, in your case the tobit regression might be overkill. That's more useful when there are both censored and uncensored values within some groups.

You, however, have two samples with exactly 0 detections out of 10 trials. Think about this as a binary detection problem. If each trial had even a 1/3 chance of leading to detection, then the probability of 0 detections out of 10 trials is only (2/3)^10 = 0.017. So you could pretty safely just say that those 2 treatments are below the detection limit and restrict a standard analysis to the remaining 3.

A 2-stage model with a binary zero/non-zero initial stage, more appropriate for count data as you note, would probably have a perfect separation problem, as two of your treatments have only 0 outcomes and all the rest have non-zero outcomes. You can exactly predict the zero/non-zero status from the category of the treatment. There are ways to deal with that, but the above suggestions for tobit regression or separating out the undetectable treatments are more appropriate here.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.