4
$\begingroup$

Here comes my case:

I conducted an experiment with roughly the following design:

30 participants, each with a unique id, were asked to rank using a likert scale how much they liked images of forests. All the participants ranked the first 8 images and then the following 5 images were randomly drawn from a pool of 15 images. Therefore, in total each participant viewed 13 images of forests, but not all images were viewed by each participant. As the responses are ordinal I have gone with a cumulative link mixed effects model to preserve the structure of the data.

And here comes when I need your expertise. So far, I believe that the random terms of my model should take into account that:

Participants (id) Image (id)

However, Im confused about how to incorporate the random effects as each participants views some but not all of the same images. Thus-far I have come to the conclusion and coded in R package ordinal it as:

clmm(likert_Rating ~ Experience + X.4 + X.3 + (1 | part_id) + (1 | Plot_ID), 
     data = TotalF)

However, no matter how many models I try, I am never sure about how to include the randoms effects using this design.

I would really appreciate if some of you could point me in the right direction as I am struggling to decide how to include the random effects.

$\endgroup$
8
$\begingroup$

Indeed, it seems that you have a crossed design, expecting that ratings from the same subjects will be correlated and that ratings for the same image will be correlated. Hence, the model you have specified seems logical. Note that the model does not require that all participants rate all images (i.e., that you have a complete balanced design). It will also work with an unbalanced design.

Check also this section of the GLMM FAQ for advice on how to code the part_id and Plot_ID variables.

| cite | improve this answer | |
$\endgroup$
5
$\begingroup$

I think that what you are doing is fine. Abstracting the situation at hand, we do not have fully crossed random effects. Nevertheless we can specify the model the same as if it were a fully crossed design. Ultimately a "deal-breaker" would be a random effect having an inadequate number of levels (usually <5) which is not the case here. One might want to check Schielzeth & Forstmeier (2008) "Conclusions beyond support: overconfident estimates in mixed models" and their subsequent 2010 paper "Cryptic multiple hypotheses testing in linear models: overestimated effect sizes and the winner's curse" for some issues that might arise by overconfident estimates and potentially misspecification of the random effects.

| cite | improve this answer | |
$\endgroup$
5
$\begingroup$

This is a fun thread with great answers! To add a bit more flesh to what was suggested here, I generated some data similar to yours using R:

part_id <- rep(1:5, each=3)
plot_id <- c(1,2,3, 1, 2, 4, 1, 2, 3, 1, 2, 5, 1, 2, 4)
rating <- c(0, 1,1, 2,1,2, 0,0,1, 2,2,1, 0, 2,1)
experience <- c(10, 12, 11, 11, 14, 13, 12, 14, 15, 11, 11, 12, 13, 15, 14)

data <- data.frame(part_id, plot_id, rating, experience)

I then converted the variables in this data set to their appropriate types:

str(data)

data$rating <- factor(data$rating, levels = c(0,1,2), labels = c(1,2,3), order = TRUE)
data$part_id <- factor(data$part_id)
data$plot_id <- factor(data$plot_id)

str(data)

Note that I recoded the original ratings (expressed as 0, 1 or 2) so they are expressed as 1, 2 or 3. Here is what the generated data look like:

   part_id plot_id rating experience
1        1       1      1         10
2        1       2      2         12
3        1       3      2         11
4        2       1      3         11
5        2       2      2         14
6        2       4      3         13
7        3       1      1         12
8        3       2      1         14
9        3       3      2         15
10       4       1      3         11
11       4       2      3         11
12       4       5      2         12
13       5       1      1         13
14       5       2      3         15
15       5       4      2         14

As you can see, there are 5 participants (whose IDs are listed in the data column named part_id). All participants see 3 images. There are 5 images in total - each participants gets to see the first two images which are labelled 1 and 2, as well as a third images selected at random out of the remaining three images which are labelled 3, 4 or 5. (The image identifiers are listed in the data column named plot_id.)

Using the xtabs() command:

xtabs(~ part_id + plot_id, data = data)

reveals the following output:

           plot_id
part_id  1 2 3 4 5
      1  1 1 1 0 0
      2  1 1 0 1 0
      3  1 1 1 0 0
      4  1 1 0 0 1
      5  1 1 0 1 0

Now, if you were to see the value 1 in all the cells corresponding to combinations of part_id and plot_id in this output, you would conclude that your two grouping factors - participant and image - are fully crossed (i.e., each participant would rate each image).
But because you see some 0's in this output, you can conclude that the two grouping factors are only partially crossed.

Here is a clmm model you can fit to the generated data:

library(ordinal)

model <- clmm(rating ~ experience + (1|part_id) + (1|plot_id), 
              data = data)

summary(model)

As pointed out in one of the other responses, the above syntax for specifying your grouping factors works even though the grouping factors are partially - not fully - crossed. (If the grouping factors were fully crossed, you would use the exact same syntax.)

The model summary reported by R is as follows:

> summary(model)
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ experience + (1 | part_id) + (1 | plot_id)
data:    data

link  threshold nobs logLik AIC   niter    max.grad cond.H 
logit flexible  15   -16.06 42.11 180(377) 1.65e-06 2.2e+04

Random effects:
Groups  Name        Variance            Std.Dev.   
plot_id (Intercept) 0.68669862403048876 0.828672809
part_id (Intercept) 0.00000000000008642 0.000000294
Number of groups:  part_id 5,  plot_id 5 

Coefficients:
           Estimate Std. Error z value Pr(>|z|)
experience  0.02879    0.39609   0.073    0.942

Threshold coefficients:
    Estimate Std. Error z value
1|2  -0.8113     4.9654  -0.163
2|3   1.1644     5.1275   0.227

The outcome value (rating) included in this model is an ordinal variable taking the values 1, 2 or 3 (as per our recoding). The predictor variable, experience, is a numeric variable taking values in the range 10 to 15.

The fitted model is really a collection of 2 sub-models (since the outcome value rating has 3 categories in total), as follows:

logit(Prob(rating for i-th participant on the j-th image <= 1)) = 
 -0.8113 - (0.02879*experience + random intercept associated with i-th participant + 
            random intercept associated with j-th image)

 logit(Prob(rating for i-th participant on the j-th image <= 2)) = 
  1.1644 - (0.02879*experience + random intercept associated with i-th participant + 
            random intercept associated with j-th image)

Note that Prob(...) is a conditional probability, since it depends on experience as well as the random intercepts included in the model.

Let's look at the predicted probabilities for a "typical" participant who rates a "typical" image (in which case we can set the random intercept associated with this participant to 0 and the random intercept associated with this image also to 0).

We know that:

logit(Prob(rating for "typical" participant on the "typical" image <= 1)) = 
   -0.8113 - (0.02879*experience) 

so that computing plogis(-0.8113 - (0.02879*experience) ) would give us Prob(rating for "typical" participant on the "typical" image <= 1), which is in effect Prob(rating for "typical" participant on the "typical" image == 1) (since our smallest rating is 1). We can plug in a value for experience in this formula, say 10, to get the desired probability:

plogis(-0.8113 - (0.02879*10))

Thus, Prob(rating for "typical" participant on the "typical" image == 1) = 0.2498898 (rounded to 0.25) when experience is 10.

We also know that:

logit(Prob(rating for "typical" participant on the "typical" image <= 2)) = 
    1.1644 - (0.02879*experience)

So we can compute Prob(rating for "typical" participant on the "typical" image = 2) as:

Prob(rating for "typical" participant on the "typical" image <= 2) - 
  Prob(rating for "typical" participant on the "typical" image <= 1)

which means that we can compute Prob(rating for "typical" participant on the "typical" image = 2) when experience is 10, say, by computing

plogis(1.1644-(0.02879*10)) - plogis(-0.8113 - (0.02879*10))

which gives us 0.4562066 (rounded to 0.46).

Because

Prob(rating for "typical" participant on the "typical" image = 3) = 
  1 - Prob(rating for "typical" participant on the "typical" image = 1) - 
      Prob(rating for "typical" participant on the "typical" image = 2)

we can compute Prob(rating for "typical" participant on the "typical" image = 3) when experience is 10 by computing

1 - plogis(-0.8113 - (0.02879*10)) - 
   (plogis(1.1644-(0.02879*10)) - plogis(-0.8113 - (0.02879*10)))

which gives us 0.2939036 (rounded to 0.29).

A convenient way to get these probabilities without computing them manually is via the ggeffects package:

 library(ggeffects)

 ggpredict(model, "experience", type="fe") 

 # Predicted probabilities of rating
 # x = experience

 # Response Level = 1

 x | Predicted |   SE |        95% CI
 -------------------------------------
10 |      0.25 | 0.23 | [-0.20, 0.70]
11 |      0.24 | 0.17 | [-0.10, 0.59]
12 |      0.24 | 0.14 | [-0.04, 0.52]
13 |      0.23 | 0.15 | [-0.05, 0.52]
14 |      0.23 | 0.18 | [-0.12, 0.58]
15 |      0.22 | 0.22 | [-0.22, 0.66]

# Response Level = 2

 x | Predicted |   SE |       95% CI
 ------------------------------------
10 |      0.46 | 0.17 | [0.12, 0.79]
11 |      0.46 | 0.17 | [0.12, 0.79]
12 |      0.45 | 0.17 | [0.13, 0.78]
13 |      0.45 | 0.16 | [0.13, 0.77]
14 |      0.45 | 0.16 | [0.14, 0.77]
15 |      0.45 | 0.16 | [0.13, 0.77]

# Response Level = 3

x | Predicted |   SE |        95% CI
-------------------------------------
10 |      0.29 | 0.28 | [-0.25, 0.83]
11 |      0.30 | 0.21 | [-0.12, 0.72]
12 |      0.31 | 0.17 | [-0.02, 0.63]
13 |      0.31 | 0.15 | [ 0.01, 0.61]
14 |      0.32 | 0.19 | [-0.04, 0.68]
15 |      0.32 | 0.25 | [-0.16, 0.81]

Adjusted for:
* part_id = 1
* plot_id = 1

Standard errors are on link-scale (untransformed).

The probabilities of interest are reported in the rows corresponding to x = 10 (where x stands for experience) and columns titled predicted for Response Level = 1, Response Level = 2 and Response Level = 3, respectively.

You can also plot the probabilities estimated by the model when experience ranges from 10 to 15 for the "typical" participant rating a "typical" image:

plot(ggpredict(model, "experience", type="fe"))

obtaining the image

enter image description here

The estimated random intercepts for participant and image can be extracted with the command:

ranef(model)

and are as follows:

> ranef(model)

$part_id (Intercept) 1 -0.35607349 2 0.56360372 3 -0.71899444 4 0.57895929 5 -0.08003278

$plot_id
              (Intercept)
1 -0.00000000000007772317
2  0.00000000000003506424
3  0.00000000000002310398
4  0.00000000000004307698
5 -0.00000000000002509980

As you can see, there is not much variation in the random effects associated with image and this is also captured in the model summary, where the standard deviation of the random intercepts associated with image is really small:

Random effects:
 Groups  Name        Variance            Std.Dev.   
plot_id (Intercept) 0.68669862403048876 0.828672809
part_id (Intercept) 0.00000000000008642 0.000000294
Number of groups:  part_id 5,  plot_id 5 

So in this case, we could probably revise the model to exclude the (1|part_id) term if need be. If we were to keep this term in the model, we could use a similar reasoning as described above to compute various probabilities for particular individuals and particular images they rated. For example,

logit(Prob(rating for participant 1 on image 1 <= 1)) = -0.8113 - (0.02879*experience + (-0.35607349) + (-0.00000000000007772317))

so that, when experience = 10,

 plogis(-0.8113 - (0.02879*10 + (-0.35607349) + (-0.00000000000007772317)))

gives us a probability of 0.3223208 (rounded to 0.32) for participant 1 on image 1 assuming their experience is equal to 10.

This tutorial might come in handy:

https://cran.r-project.org/web/packages/ordinal/vignettes/clmm2_tutorial.pdf.

in terms of understanding what is being modelled in similar models (but which include a single random intercept).

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ This is really a great walk through of these ideas. Thank you. $\endgroup$ – MrDrFenner Feb 6 at 16:26
  • $\begingroup$ You're welcome, @MrDrFenner! I was beginning to think that going into such detail was not warranted given the lack of reaction to the answer. 😛 $\endgroup$ – Isabella Ghement Feb 6 at 17:07
  • 1
    $\begingroup$ It may have been overkill ... but I might encourage you to post this as a blog somewhere .... or ask-and-answer a question here on StackExchange ... for which this great answer will then have a perfect question. $\endgroup$ – MrDrFenner Feb 6 at 17:28
  • $\begingroup$ I noticed that my (rare) blog posts receive little traction - there is not much appetite for detailed posts out there and I really don't have the time to write them only to see them languish without much impact. So this will have to do. $\endgroup$ – Isabella Ghement Feb 6 at 18:08

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.