Animal model repeated measures in brms: levels of within group covariance matrix do not match names of grouping levels

Interfaces brms
, ,
1
  • Operating System: Microsoft Windows 11
  • brms Version: 2.16.1

Hello!

I am having the same problem that Rhainer had, here. It was never resolved, and so I am writing to the brms community for help.

I’m trying to fit a repeated measures animal model to estimate repeatability. In my simple model, I’d include a behavioral response variable Resp, fixed factor Sex, relatedness matrix covmat, and individual ID. Normally, when specifying brms models, I’d link ID to the animal covmat using the ‘gr’ term:

model <- brm(
  Resp ~ Sex + (1|gr(ID,cov=A)),
  data=data,
  data2=list(A=covmat),
  family = gaussian(),
  warmup = 100, iter = 5000,

However, when trying to incorporate a second ID random effect term, to account for fixed among-individual differences (permanent environment effects), this is where my models fail. Below is my very simple model incorporating ID_pe and the error message I receive:

model <- brm(
  Resp ~ Sex + (1|gr(ID,cov=A)) + (1|ID_pe), # ID_pe is the permanent environmental effects
  data=data,
  data2=list(A=covmat),
  family = gaussian(),
  warmup = 100, iter = 5000,
  chains = 1, cores = 1,thin=1)


Error: Levels of the within-group covariance matrix for 'ID' do not match names of the grouping levels.

As Rhainer did in his post, I checked to ensure my group names matched between my covariance matrix and my phenotype data–and they do.

One possible solution is unlinking the animal covmat from the animal ID term. When I run my model as specified in this tutorial, the model runs. Code:

model <- brm(
  Resp ~ Sex + (1|ID) + (1|ID_pe), # no grouping term linking ID to animal covmat
  data=data,
  data2=list(A=covmat),
  family = gaussian(),
  warmup = 100, iter = 5000,
  chains = 1, cores = 1,thin=1)

My questions are therefore:

  1. What could be causing the mismatch of group level names?

And perhaps I don’t need to know the answer to number 1 if:

  1. Is this second way of model specification equivalent to the first? That is, by not explicitly linking ID to the covariance matrix covmat, will I still be able to estimate additive genetic variance?

Much thanks to the community in advance,
Tracy

2

Howdy!
I do not know the answer to your questions (I don’t work with phylogenetic models), but it seems that the post that you link to was not resolved in part because no reproducible example was provided. Perhaps you could provide the data, or simulated data, and the model code, so that someone could reproduce your problem.

Also, for question 2, you can use make_stancode() for your brms model and look at the Stan code directly, to see if the second specification encodes the model and assumptions that you are looking for.

1 Like
3

Question 1:
Hi. I have been fitting many models like these, i.e. of the form trait ~ 1 + (1 | animal) + (1 | gr(animal_A, cov = A)) without issues.

Also note that this is basically equivalent to the model with measurement error in the brms phylogenetics vignette, which works fine.

So I’m not sure what is going wrong, a reproducible example would help!

Question 2:

No, I don’t think these are the same models. I think A is just ignored. We can borrow some code from the phylogenetics vignette:

library(brms)

phylo <- ape::read.nexus("https://paul-buerkner.github.io/data/phylo.nex")
data_repeat <- read.table(
  "https://paul-buerkner.github.io/data/data_repeat.txt",
  header = TRUE
)
data_repeat$spec_mean_cf <- with(data_repeat, sapply(split(cofactor, phylo), mean)[phylo])
A <- ape::vcv.phylo(phylo)

model_repeat1 <- brm(
  phen ~ spec_mean_cf + (1|gr(phylo, cov = A)) + (1|species),
  data = data_repeat,
  data2 = list(A = A),
  chains = 2, cores = 2,
  iter = 2000, warmup = 1000
)

model_repeat2 <- brm(
  phen ~ spec_mean_cf + (1|phylo) + (1|species),
  data = data_repeat,
  data2 = list(A = A),
  chains = 2, cores = 2,
  iter = 2000, warmup = 1000
)

As one might expect, the second model struggles, since it cannot differentiate the phylo from species effects, and the estimates are different between the models.

2 Likes
4

@jd_c and @Ax3man,

I can confirm that these aren’t the same models–just got feedback from the authors of the tutorial, and their code was a mistake that they’ve now fixed (that is, you do indeed need to specify the “gr” term).

Yes–seems like a reproducible example would be helpful here! I’m attaching some simulated data with animal IDs. It seems like my IDs for my covariance matrix and repeated measures phenotype data must not match…but I can’t seem to figure out how to rectify this, and I don’t have this problem when specifying a simple univariate model with averaged data.

Could it be that not all individuals have repeated measures? (That is, some individuals only have 1 measure, and they somehow get “dropped”, but the covariance matrix doesn’t drop them and thus there is a mismatch??)

Thanks for your replies!
genetic covariance matrix (from SNPs) sim_covmat.csv (480.9 KB)
simulated phenotype data sim_data.csv (21.5 KB)

5

Ok, I can now reproduce, as follows:

library(brms)

A <- read.csv('~/Downloads/sim_covmat.csv')
# first column contains the rownames, so:
A <- tibble::column_to_rownames(A, 'X')
A <- as.matrix(A)

d <- read.csv('~/Downloads/sim_data.csv')

m <- brm(
  Resp ~ 1 + (1 | gr(animal, cov = A)),
  data = d, data2 = list(A = A),
  chains = 2, cores = 2, iter = 1000
)

Gives

Error: Levels of the within-group covariance matrix for 'animal' do not match names of the grouping levels.

However, checking whether it is true that all animals are in A, that is not the case:

setdiff(d$animal, rownames(A))
# [1] "SGTf03"  "SGTf10"  "SGTf15"  "SGTf27"  "SGTm062" "SGTm080" "SGTm082"
# so those are in d, but missing from A

'SGTf03' %in% d$animal
# TRUE
'SGTf03' %in% rownames(A)
# FALSE
'SGTf03' %in% colnames(A)
# FALSE

Dropping those individuals from the data, the model works:

d_reduced <- d[d$animal %in% rownames(A), ]

m <- brm(
  Resp ~ 1 + (1 | ID) + (1 | gr(animal, cov = A)),
  data = d_reduced, data2 = list(A = A),
  chains = 2, cores = 2, iter = 1000
)

So brms appears to be behaving correctly here. But leaves the question, why are your A matrices missing individuals?

No, individuals with one observation don’t get dropped.

4 Likes
6

Thanks so much for checking on my problem!! What a simple solution that I missed. I’d insert a “facepalm” emoji here if they were readily available…

So brms appears to be behaving correctly here. But leaves the question, why are your A matrices missing individuals?

The answer is simply that I measured data from animals in the wild, and not all animals had both genotype and phenotype recorded! In my previous analyses, I had used univariate models on averaged phenotypes, and I forgot that in the averaged-phenotype dataset, I had only included animals that also had genotype data, which is why they worked.

Thanks so much for checking this out for me and for being unjudgmental about it. I’ve never actually used the “setdiff” command in R, and I will remember this one for the future.

1 Like
7

Great, yes, that makes a lot of sense.

And no worries, good luck with your analyses!

1 Like