I have the following situation: General practitionar (gp), patient (pat) and consultation (cons).

Each gp has several patient and each patient can have 1 and more consultations with a specific characteristics cons_x. The outcome is y=0/1 on level consultation.

A possible model could be (in lmer/R notation)

y ~ gp_sex + pat_sex + pat_age + cons_x + (1|gp) + (1|gp:pat)

gp is a random variable and pat is nested in gp. Since the outcome is 0/1 (referral to specialist yes/no) the suitable model is a multilevel logistic regression.

(1) The problem: 50% of patients have only one consultation, that is, there is no repetition, others have 2, 3 until 6 consultations.I would like to know how to cope with this kind of situation. With other words: I have groups (patients) with only one consultation. So, within the group is perfect correlation due to the only member in the group. Does the mixed models framework cope with this?

(2) Another question which arises: The number of consultations weights the effect on patient since the same patient has as many entries as he has consultations. The more consultation a patient has the more patient characteristics are considered. This wouldn't be the case if all patient would have the same number of consultations.

Every help and hint is appreciated.

To simplify the analysis I could aggregate on level patient: If a patient has several consultations and at least one consultation has 1 then the patient has outcome 1 otherwise 0. Then I have a simpler model:

y ~ gp_sex + pat_sex + pat_age + (1|gp)

In this case problem (1) will be obsolete but than still I have problem (2) since each gp has different numbers of patients.

  • 1
    $\begingroup$ The short answer is that you basically don't need to worry. The weighting gets taken care of automatically. If there are a very small number of patients with multiple consultations then you might have a numerical problem, but otherwise everything should be fine. $\endgroup$ – Ben Bolker May 12 '14 at 22:44
  • $\begingroup$ @Ben. Thank you very much for your comment. That's what I hoped to hear. I would like to give you the bounty but obviously the term is exceeded. $\endgroup$ – giordano May 14 '14 at 21:17

It's tough without truly understanding what your outcome y is, but I'll give it my best shot.

(1) This sounds like it might be approximated by a λ=1 Poisson distribution. Unless you're seeing it messing up the variance components of the model, I wouldn't worry about it. If it is, I would try transforming it by sqrt(cons_x) link and seeing if that works better. Don't forget to check your model for sensibility by plugging in some examples in your final model, varying the # of consultations, and seeing if it makes sense.

(2) Sounds like interaction variables need to be added to see the actual effect of consultations on the sex and age.

y ~ gp_sex + pat_sex + pat_age + cons_x + cons_x*pat_sex + cons_x*pat_age + cons_x*pat_age*pat_sex + ....
  • $\begingroup$ Thanks for answer. The outcome is dichotomous (0/1) as it is described in the text. It can be, for example, referral to a specialist. Sorry that I didn't make it clearer. With this outcome the model to use is a logistic regression. But whatever the outcome is the two questions are not yet answered. I will update my questions to make them clearer. $\endgroup$ – giordano May 11 '14 at 15:50
  • 1
    $\begingroup$ Your clarification really helps, thank you. It's beginning to sound more like a data structure issue. I agree this is a logistic regression now that you state the outcome as "refer to a specialist". Your variables are as follows: gp, pat, pat_age, pat_sex, cons_x. Your model is y ~ gp:pat + pat_age + pat_sex + cons_x + cons_x*pat_sex + cons_x*pat_age + cons_x*pat_age*pat_sex. You are going to need to collapse your multiple records due to consultations into one record, where the variable cons_x indicates the # of consultations (use R to count the # of records per patient). $\endgroup$ – Gary Chung May 12 '14 at 21:18
  • $\begingroup$ Thanks for response. In the update of my question I collapsed to patient and the consultation characteristics disappears. It's a good idea to use the number of consultation as variable. What I do not anderstand is gp:pat: in R terms it means interactions. If this is what you mean I do not understand the model. Anyway, obviously I don't have to care as Ben Bolker states (see comment). $\endgroup$ – giordano May 14 '14 at 21:07
  • 1
    $\begingroup$ Ben's answer is correct but it is answering (1). In (2), you state "The more consultation a patient has the more patient characteristics are considered." I take this to mean that pt characteristics affect the outcome differently when the # of consultations vary. This is an interaction between pat_age pat_sex and cons_x. My apologies on the confusion with gp:pat, I wanted that to mean pat nested within gp, not interaction. The cons_x*pat_sex and so forth is intended to mean an interaction effect. $\endgroup$ – Gary Chung May 15 '14 at 19:45
  • $\begingroup$ Thanks for help. I appreciated your help and advices. Next time I will try to make clearer statements. $\endgroup$ – giordano May 3 '16 at 6:16

I'm certainly no expert and would love others to comment on this, but:

I'm not sure what your outcome is but you said it was measured 1/0. for the sake of the example i'll pretend that your outcome is "happy with consultation" yes=1 and No=0.

I think that for the people with multiple consultations you have a easier task. You could work out proportion happy responses for each person based on their multiple consultation (e.g person 1 had 6 consultations was happy 3 times and not happy 3 time = 50% happy). You could then work out an 'Overall'average happiness across all people by adding up the % happy for each person and dividing by the number of people. You would then get an 'overall average % happy' measure. You'd also get a standard error around this overall % . You basically now have data that could be analysed with a t test or other measures based on a normal distribution.

Where you don't have a repeat for the person, you cannot do this. the best you can do is work out a proportion of happy responses, across all participants. lets say you had 4 participants who only had 1 consultation (p1 was happy=1 p2 was happy=1, p3= 1 and p4=0) so the proportion happy is 75%. You would now have to use binomial tests e.g a z test - which i believe are not as powerful (I have just posted a question on this).

YOu may have to analyse your data for the no repeats, and the repeats separately. Unfortunately you can't just pretend that you repeated data is from different people because this would violate assumption of independence, so i sense that a different approach is needed with the repeated measures vs no repeated measures data.

  • $\begingroup$ Thanks for your detailled comment. I understand your intention and this kind of analysis was already performed but the disadvantage is that the hierarchical structure get lost. I would also like to measure the contribution of the gp. I thought that mixed models can cope within certain limits with single observations since there are still pat with several observations. $\endgroup$ – giordano May 8 '14 at 9:59

I can deliver an answer to my question only empirically using a simulation. Using this two cross validation contributions mixed and logistic I could create some fake datasets and using the mixed logistic regression. With R and glmer from library(lme4) I used this formula:

 fit1 <- glmer(y ~ x1 + (1|j), data = d, family=binomial)

y is a dichotomous variable, x1 is continuous and j is random. First I built a balanced dataset d with the grouping variable j with 20 groups. Then I build two datasets from d with m groups which have only one observation.
The simulations shows that the variance of j and the fixed coefficient b1 of x are very similar to the "true" values for all kind of datasets. This is due to "perfect" randomization. In reality there will be some bias.
The limits of this answer: lacking of theoretical foundation (but intuitively it makes sense). The simulation could be improved if packed in a loop to get many estimates and comparing the average with the true values var(uj) and coefficient b1.


# -- Model
# yj[i] = b*0 b1*xj[i]
# b0 = g00 + u0j, u0j ~ N(0,1)
# b1 = const
# => zj[i] = g00 + u0j[i] + b1*xj[i]

# -- Libraries

# -- Create balanced dataset d
# Number of clusters (level 2)
N <- 20
# Number of observations (level 1) for cluster j
nj <- 200
# intercept
g00 <- 1
# slope
b1 <- 3
# Vector of clusters indices 1,1...n1,2,2,2,....n2,...N,N,....nN
j <- c(sapply(1:N, function(x) rep(x, nj)))
# Vector of random variable
uj <- c(sapply(1:N, function(x)rep(rnorm(1), nj)))
# Vector of fixed variable
x1 <- rep(rnorm(nj),N)
# linear combination
z <- g00 + uj + b1*x1
# pass trhough an inv-logit function
pr <- 1/(1 + exp(-z))
# bernoully response variable
y <- rbinom(N*nj,1,pr)
d <- data.frame(j, y=y, z=z,x1=x1, uj=uj)

# -- Create unbalanced datasets d2 and d3
# sort j
d <- sqldf("SELECT * FROM d ORDER BY j")

# count each observation within j
d$ord    <-  NA
d$ord[1] <- 1
k <- 2
for (i in 2:nrow(d) ) {
    if ( d$j[i] == d$j[i-1] ) {
      d$ord[i] <- k
      k = k+1
  } else {
      d$ord[i] <- 1
      k = 2
# result

# Define a sample with m groups which have only one observation
m <- 5
d2 <- subset(d, j %in% c(1:m)
                | ( j %in% c((m+1):20) & ord == 1)
# Another sample
d3 <- subset(d, j %in% c(1:m)
                | ( j %in% c((m+1):20) & ord == 10)

# Fit with balanced dataset d

# fit with unbalanced data d2
fit2 <- glmer(y ~ x1 + (1|j), data = d2, family=binomial)

# fit with unbalanced data d3
fit3 <- glmer(y ~ x1 + (1|j), data = d3, family=binomial)

Your Answer

