1
$\begingroup$

I have a dataset with one continuous response variable (time), a 'treatment' explanatory variable and 5 other fixed factors:

summary(dat)
area         date    day_time    treatment       time        trials  
A:31   05.04.14: 8   A:23     Control :32   Min.   :0.0000   Orca2_1 :  2  
B: 8   06.04.14: 8   B:28     Orca2 ON:15   1st Qu.:0.0000   Orca2_14: 2  
C:26   07.04.14: 8   C:14     SR2 ON  :18   Median :0.0278   Orca2_16: 2  
       08.04.14: 8                          Mean   :0.1272   Orca2_17: 2  
       25.04.15: 6                          3rd Qu.:0.2023   Orca2_18: 2  
       28.04.15: 6                          Max.   :0.9216   Orca2_19: 2  
       (Other) :21                                           (Other) :53  

I would like to build a linear mixed effect model like:

M1 = lmer(time ~ treatment + (1|date) + (1|area) + (1|trials) + (1|day_time) , data=dat)

However, the data is zero-truncated as seen on the histogram: histogram of time

When I run the model, I get an abnormal distribution of residuals:

norma qqplot

I tried to transformed the time variable with an inverse hyperbolic sine transformation:

transformed_time <- log(dat$time + sqrt(dat$time^dat$time +1))

and then run the model again:

M2 = lmer(transformed_time ~ treatment + (1|date) + (1|area) + (1|trials) + (1|day_time) , data=dat)
qqnorm(residuals(M2))

It gets a bit better but still not convincing.

I am getting confused with everything that I am reading on zero-truncated models for continuous data... I was wondering if someone could help me find a way to build a valid model.

Here is the data:

EDIT : added a time_cat column for time categories A-I

> dput(dat)
structure(list(area = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 1L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
3L, 3L, 3L, 3L, 1L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("A", "B", 
"C"), class = "factor"), date = structure(c(1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 6L, 6L, 
6L, 8L, 7L, 9L, 10L, 10L, 10L, 11L, 12L, 12L, 12L, 3L, 3L, 3L, 
4L, 4L, 4L, 4L, 8L, 10L, 10L, 10L, 11L, 12L, 12L, 12L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 9L, 9L, 9L
), .Label = c("05.04.14", "06.04.14", "07.04.14", "08.04.14", 
"11.04.14", "25.04.15", "26.04.15", "26.05.15", "27.04.15", "28.04.15", 
"29.04.15", "30.04.15"), class = "factor"), day_time = structure(c(1L, 
1L, 2L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 1L, 2L, 2L, 3L, 
1L, 2L, 1L, 2L, 3L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 
2L, 3L, 1L, 2L, 2L, 3L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 
2L, 3L, 1L, 2L, 3L, 3L, 1L, 2L, 2L, 1L, 2L, 3L, 1L, 1L, 1L, 2L
), .Label = c("A", "B", "C"), class = "factor"), treatment = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L
), .Label = c("Control", "Orca2 ON", "SR2 ON"), class = "factor"), 
    time = c(0.2148, 0.1814, 0.106300000000005, 0.248799999999999, 
    0.129899999999999, 0.109099999999998, 0, 0.145200000000003, 
    0.1522, 0, 0.202300000000001, 0.921599999999998, 0.580300000000001, 
    0.1327, 0.617799999999995, 0.3309, 0.3127, 0.311299999999999, 
    0.151499999999999, 0, 0, 0, 0, 0.0806000000000004, 0.262699999999995, 
    0, 0, 0, 0, 0.2224, 0, 0, 0.136900000000004, 0.0743999999999971, 
    0, 0, 0.0784999999999982, 0.360700000000001, 0, 0, 0.0277999999999992, 
    0, 0, 0, 0, 0, 0, 0.238399999999999, 0.169600000000003, 0.394000000000005, 
    0, 0, 0, 0, 0.152200000000008, 0.151499999999999, 0.440600000000003, 
    0.331499999999998, 0, 0, 0, 0, 0, 0.296800000000005, 0), 
    time_cat = structure(c(4L, 3L, 3L, 4L, 3L, 3L, 1L, 3L, 3L, 
    1L, 4L, 9L, 7L, 3L, 8L, 5L, 5L, 5L, 3L, 1L, 1L, 1L, 1L, 2L, 
    4L, 1L, 1L, 1L, 1L, 4L, 1L, 1L, 3L, 2L, 1L, 1L, 2L, 5L, 1L, 
    1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 4L, 3L, 5L, 1L, 1L, 1L, 1L, 
    3L, 3L, 6L, 5L, 1L, 1L, 1L, 1L, 1L, 4L, 1L), .Label = c("A", 
    "B", "C", "D", "E", "F", "G", "H", "J"), class = "factor"), 
    trials = structure(c(16L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 
    33L, 33L, 1L, 7L, 11L, 12L, 13L, 14L, 15L, 17L, 19L, 20L, 
    21L, 22L, 23L, 2L, 25L, 3L, 4L, 5L, 6L, 8L, 9L, 10L, 1L, 
    7L, 11L, 12L, 13L, 14L, 15L, 2L, 3L, 4L, 5L, 6L, 8L, 9L, 
    10L, 16L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 17L, 18L, 19L, 
    20L, 21L, 22L, 23L, 24L, 24L, 25L), .Label = c("Orca2_1", 
    "Orca2_14", "Orca2_16", "Orca2_17", "Orca2_18", "Orca2_19", 
    "Orca2_2", "Orca2_20", "Orca2_21", "Orca2_22", "Orca2_3", 
    "Orca2_4", "Orca2_5", "Orca2_6", "Orca2_7", "SR2_1", "SR2_10", 
    "SR2_11", "SR2_12", "SR2_14", "SR2_15", "SR2_16", "SR2_17", 
    "SR2_18", "SR2_19", "SR2_2", "SR2_3", "SR2_4", "SR2_5", "SR2_6", 
    "SR2_7", "SR2_8", "SR2_9"), class = "factor")), .Names = c("area", 
"date", "day_time", "treatment", "time", "time_cat", "trials"
), class = "data.frame", row.names = c(NA, -65L))
$\endgroup$
1
$\begingroup$

You could try glmer, which allows a variety of non-normal distributions.

Your data are not "zero-truncated" -- your response is time which is perforce non-negative. No truncation happened. It's just not normally distributed data.

Do the zero time observations correspond to genuinely instantaneous events? Or do you actually have grouped observations .. whereby t=0 actually means 0 < t < short time interval?

Later: the OP has supplied raw data, so in light of that, I have changed my answer.

The data set is not very large (65 observations), divided almost equally between 0's and non-zeros (32, 33). The non-zero observations are very long-tailed. I'm not sure how much analysis a data set like this can support. As a first crack, you can dichotomize time to 0 or 1, according as time is 0 or greater than 0. Then build up a model using logistic regression.

Note that most of these variables are not random effects. Treatment, day_time and Area should be fixed. Trial is possibly random, but possibly also irrelevant.

I ran the following:

dat$time2 <- ifelse(dat$time ==0, 0, 1)
M1 <- glm(time2 ~ area + day_time + treatment, data=dat, family=binomial())
summary(M1)

M2 <- glmer(time2 ~ area + day_time + treatment + (1|trials), data=dat, family=binomial())
summary(M2)

In model M2, the variance of random effect trials was 0, so there doesn't seem to be much difference between trials as such.

Here are the results of model M1.

Call:
glm(formula = time2 ~ area + day_time + treatment, family = binomial(), 
    data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9621  -0.6714   0.3580   0.6508   2.5680  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)         1.7673     0.7667   2.305   0.0212 *  
areaB              -1.1782     1.0218  -1.153   0.2489    
areaC              -3.1522     0.7882  -3.999 6.36e-05 ***
day_timeB           0.9479     0.7688   1.233   0.2176    
day_timeC          -0.3229     0.9400  -0.343   0.7312    
treatmentOrca2 ON  -1.8748     0.8734  -2.147   0.0318 *  
treatmentSR2 ON    -0.9380     0.7753  -1.210   0.2264    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 90.094  on 64  degrees of freedom
Residual deviance: 62.176  on 58  degrees of freedom
AIC: 76.176

Number of Fisher Scoring iterations: 5

It looks as if Area C matters and treatment Orca2.

Note that glm() is part of base R and glmer() is in package lme4.

I'm not sure it's worth doing much more than this on your data. If you want to go further, you need to treat this as a ZIP design. The non-zero stuff seems to be a gamma, or something like it. As a first pass, you could try and model a gamma to the non-zero part of the data and see if what you get makes sense.

In Poisson ZIP's, the zeroes belong either to the Poisson, or to the "inflation". When the non-zero part is a continuous distribution, strictly speaking, you can't model a gomp of zeroes, so you would need to jitter the data away from zero. I guess the question you need to ask yourself is what you believe the treatments are going to do to the times? Do they basically change the number of zeroes (in which case my binomial kludge is good enough), or will they change the non-zero times? In any case, I don't think you have enough data to go too deeply into this. If you had more, I would suggest fitting two models: a binomial model that looks at the zero/non-zero dichotomy, and something else (possibly gamma) to look at the non-zero values on their own.

$\endgroup$
  • $\begingroup$ Thanks a lot for your advice Placidia. The zeros in this data set are genuine zeros. I will try a glm with time intervals, that sounds like a good idea! I will let you know the outcome. $\endgroup$ – user3406207 Aug 19 '15 at 0:55
  • 1
    $\begingroup$ So I have now a new variable 'time_cat' with 9 factors A-I. However I am not too sure how to put that in a nb glm as it is not a count variable yet. I don't want to sum it up either as I will loose the information from the other variables... Can you help? (ps I edited the dataset in my question) $\endgroup$ – user3406207 Aug 19 '15 at 4:23
  • $\begingroup$ I looked at your data. Negative binomial was an error on my part. I changed my answer. Hope this helps. $\endgroup$ – Placidia Aug 19 '15 at 17:40
  • $\begingroup$ Brilliant. This all makes sense to me (and to the data). As you said, the treatments should only change the number of zeros and not influence the non-zero too much. The binomial model is enough for what I need and I accept it gladly! Thank you so much for the help. $\endgroup$ – user3406207 Aug 20 '15 at 2:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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