Repeated measures with single measurements

Question

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.

Update:
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.

Answer 1

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 + ....

Answer 2

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.

Answer 3

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.

Simulation:

# -- 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
library(lme4)
library(sqldf)

# -- 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
d[c(190:210),]

# 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
fit1 <- glmer(y ~ x1 + (1|j), data = d, family=binomial)
summary(fit1)

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

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