I have data that look like the example below. There are 3 different groups (g1, g2 and g3) of subjects which needed either 1, 2 or 3 procedures (p1, p2, p3) respectively according to the course of the disease. Hence, there is only p1 for the group g1; p1, p2 for the group g2; p1, p2, p3 for the group g3. There are several measurements (8 in total for the real data): some of which are continuous (such y1 in the example below) other dichotomous (such as y2 in the example below), with missing data (NA). Each measurement was performed at the time of each procedure. Moreover there is the date at which each procedure was performed and the timespan in days between each consecutive procedure of each procedure. I'm trying to compare the mean of y1 and the frequency of y2 between the several time points of the procedure.

Here is a representative simulation of the data:


# define the number of subjects in each groups
n_g3 = 20
n_g2 = n_g3 * 2
n_g1 = n_g3 * 5

# 3 different groups
id_p1 = paste0("ID",1:c(n_g1 + n_g2 + n_g3))
id_p2 = paste0("ID",c(n_g1+1):c(n_g1 + n_g2 + n_g3))
id_p3 = paste0("ID",c(n_g1 + n_g2+1):c(n_g1 + n_g2 + n_g3))
id = append(append(id_p1, id_p2), id_p3)

# 3 different groups
groups = c(rep("g1", n_g1), rep("g2", n_g2), rep("g3", n_g3), rep("g2", n_g2), rep("g3", n_g3), rep("g3", n_g3))

# 1st, 2nd or 3rd procedure
procedure = c(rep("p1", n_g1), rep("p1", n_g2), rep("p1", n_g3), rep("p2", n_g2), rep("p2", n_g3), rep("p3", n_g3))

# date of procedure and timespan between procedures
date_p1 = as.Date(as.Date(rep("2012-01-30",c(n_g1 + n_g2 + n_g3))) - replicate(c(n_g1 + n_g2 + n_g3), 40*rnorm(1)) %>% round())
timespan_p1_to_px = c(rep(0, n_g1 + n_g2 + n_g3), 
                  replicate(c(n_g2),500 +(40*rnorm(1))) %>% round(digit=0),
                  replicate(c(n_g3),650 +(40*rnorm(1))) %>% round(digit=0),
                  replicate(c(n_g3),1500 +(40*rnorm(1))) %>% round(digit=0))

date_construct = c(date_p1, date_p1[(n_g1+1):(n_g1 + n_g2 + n_g3)], date_p1[(n_g1 + n_g2+1):(n_g1 + n_g2 + n_g3)])

date_px = as.Date(as.Date(as.Date(date_construct)) + timespan_p1_to_px)

# measurement n°1 (continuous)
y1 = c(
  (replicate(c(n_g1 + n_g2 + n_g3), 20+(3*rnorm(1))) %>% round(digit = 1)) %>% abs(),
  (replicate(c(n_g2 + n_g3), 25+(3*rnorm(1))) %>% round(digit = 1)) %>% abs(),
  (replicate(n_g3, 26+(3*rnorm(1))) %>% round(digit = 1)) %>% abs()
    y1[rbinom(c(n_g1 + n_g2 + n_g3),1,0.2) == 1] = NA

# measurement n°2 dichotomous
y2 = c(
  (rbinom(c(n_g1 + n_g2 + n_g3),1,0.8)),
  (rbinom(c(n_g2 + n_g3),1,0.5)),

y2[rbinom(c(n_g1 + n_g2 + n_g3),1,0.2) == 1] = NA

data = data.frame(id, groups, procedure, date_px, timespan_p1_to_px, y1, y2)

### create subset of data #####

# subset p1 to make the comparisons of all the variable at baseline
data_p1 = data[which(data$procedure == "p1"),]

# subset g2 to make comparisons within the group with 2 procedures
data_g2 = data[which(data$groups == "g2"),]

# subset g3 to make comparisons within the group with 3 procedures
data_g3 = data[which(data$groups == "g3"),]


I have the 5 following questions:

Question N°1: I would like to compare the baseline measurement at p1 for g1, g2 and g3. Those are three different groups with all different individuals. I would compare the mean of each group with an ANOVA as follow:

Question N°2: I would like to compare if there is a difference in the measurements within the group g2 between time p1 and p2. For this I think I should use a linear mixed model. I would use the following command:

Question N°3: I would like to compare if there is a difference in the measurements within the group g3 between time p1 and p2 and between p2 and p3. For this I think I should use a linear mixed model. I would use the following command:

Question N°4: How should i deal with a dichotomous variable of time y2. should I use the following lines:

Question N°5: Finally: Would it be possible to accound for the timespan elapsed between p1 and p2 in g2 and g3, and between p2 and p3 in the group g3? if yes what would be the line of code

  • 1
    $\begingroup$ Using a mixed model would indeed be appropriate in this case because it would give you unbiased results under the missing at random missing data mechanism (provided that the random-effects structure is correctly/flexibly specified). From your description it is was clear if you want to separately analyze y1 and y2 or jointly (e.g., are you interested in the association between y1 & y2?). And moreover why you have missing data (i.e., is it missing completely at random or missing at random)? $\endgroup$ – Dimitris Rizopoulos Aug 31 '18 at 10:55
  • $\begingroup$ Thank you @DimitrisRizopoulos for your answer. So to answer your questions: 1. I am not particularly interested in association between y1 and y2. What I want to do is if there is relevant difference between y1_t1 and y1_t2 and y1_t3 (as well as for y2) between the different time points. 2. I have missing data that I assume being completely at random. Do you know some code (e.g. using lme4) that could model this question? $\endgroup$ – ecjb Aug 31 '18 at 14:26
  • $\begingroup$ @DimitrisRizopoulos. I changed the display of the data and made the questions more precise with an attempt to answer them $\endgroup$ – ecjb Sep 5 '18 at 20:12

Answer to question N°1:

Using an ANOVA:

mod1 = aov(formula = y1~groups, data = data_p1)

or using a linear model and the multcomp package as described in this stackoverflow question:

mod1_2 = lm(y1~groups, data = data_p1)
mod1_2 = glht(mod1_2, linfct = mcp(groups = "Tukey"))

Answer to question N°2:

Using the lmerTest to get the p value in summary of the regression:

mod2 = lmer(y1 ~ procedure + (1|id), data = data_g2)

Answer to question N°3:

Performing a linear mixed models with id as random variable and performing a Tukey post-hoc test to look at all the between groups comparison using the multcomp package, as shown on this stack exchange question

mod3 = lmer(y1 ~ procedure + (1|id), data = data_g3)
summary(glht(mod3, linfct = mcp(procedure = "Tukey")), test = adjusted("holm"))

Answer to question N°4:

mod4 = glmer(y2 ~ procedure + (1|id), data = data_g2, family = binomial)

mod4_2 = glmer(y2 ~ procedure + (1|id), data = data_g3, family = binomial)

Answer to question N°5:

I'm really not sure about this one

mod5 = lmer(y1 ~ timespan_p1_to_px + (1|id) + (1|procedure), data = data_g3)

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