Could BRMS be used to estimate an ellipse? (Issues with estimating unseen parameters)

I have collected some experimental data in which participants focus on a centre and have to identify surrounding objects. It has long been identified in the literature that some positions may be easier to identify than others. I have therefore made a staircase procedure that tries to move the objects back and forth until the accuracy is roughly equal amongst positions:

|309px;x309px;

I have already used BRMS to compare models with different assumptions (for instance an upper/lower division vs a left/right division). However, I have not had success with modelling an underlying ellipsoid distribution. Investigating wiki, it seems like one relevant equation could be:

image
I believe it could be possible to estimate with rstan but I wonder if it would be possible to estimate using BRMS. I would much prefer to stay in this environment as I have already done most of the analysis here and it would be great to formally compare these models with loo to see which would be best to describe the data.

I have created a script which simulates ellipsoid data and attempts to analyse it here.

My thought of how a model should be specified should be something like:

brm(data = dat,
family = gaussian,
mvbind(xposa(cos(t)),
yposb(sin(t))) ~ 1 +(1|part), chains = 2, cores = 2, iter=500)

Which returns:

Error: The following variables are missing in ‘data’:
‘a’, ‘t’, ‘b’

I have seen unknown variables in rstan before but not in BRMS, is this possible?

My best (very sorry attempt) does not estimate the data very well:

brm(data = dat,
family = gaussian,
mvbind(xpos, ypos) ~ 1 +(1|part), chains = 2, cores = 2, iter=500)

Is there a way that I can progress without having to start the analysis from scratch? All your response will be appreciated, I am a bit stuck here.

1 Like

Hi,
sorry for taking to long to respond. I am however not sure I follow the question exactly. What do you mean by “ellipsoid distribution”? Are you trying to determine the parameters of the ellipse that would best fit your data? Because in your code to generate data, you simply add normal noise to the equation of ellipse - does the simulated code match the data generating process you imagine?

However, if the experiment actually involved only moving points closer to and farther away from the centre, then I would not work with the x and y positions, but rather with the distances from the center. So you would have the mean distance d as a function of angle t:

d(t) = \sqrt{a^2 \cos^2(t) + b^2 \sin^2(t)}

Now d(t) would be the mean of a normal distribution around which you expect your observations to lie, e.g.:

y_i \sim N(d(t_i))

EDIT: Since the distance is always positive, using log normal distribution (family = "lognormal") is likely preferrable.

With some struggle, this should be possible to achieve using the non-linear modelling features in brms. Something along the lines of (code not tested, just a sketch):

  bf(distance ~ sqrt(a ^ 2 * sin(angle)^2 + b^2 * cos(angle)^2),
     a ~ 1, b ~ 1, 
     nl = TRUE)

where distance and angle are data. This would estimate a single ellipse for all data. You could then have a and b depend on part to add the hierarchical structure. Also note that non-linear formulae behave a bit differently than the “standard” brms formulae, so read the vignette carefully. This parametrization also makes it visible that there is a symmetry in the ellipse, so you would need to constraint a and b to be positive.

Hope any of this is making sense.

2 Likes

Thank you for your response, Martin. Please don’t be sorry about the response time, I am happy to wait for such helpful input. It would maybe be easier for the both of us, if I make the experimental data available. Here is a minimal data set containing:
Distance, Angle & Part (unique participant)

While I did get a non-linear model up and running from your advice it seems that something is still off. I specify the model as:

priors <- c(set_prior(“normal(0,2)”, nlpar = “a”, lb = 0),
set_prior(“normal(0,2)”, nlpar = “b”, lb = 0))

nlModel <- brm(bf(Distance ~ sqrt(a ^ 2 * sin(Angle)^2 + b^2 * cos(Angle)^2),
a ~ 1, b ~ 1,
nl = TRUE), data = data, prior = priors, cores = 2, chains = 2)

And that returns the following estimate:

The code that has been analysing it thus far in a linear model looks like this:

lModel <- brm(Distance ~-1+anglesF +(1|Part), data = data, chains = 2, iter=2000, cores = 2)

And this is the estimated shape:

It thus seems that the estimated non-linear model fits the data quite poorly and I guess that the shape is strictly speaking not an ellipse. A comparison with loo also has strong preference for the linear model (Though note that I haven´t added random effects for participants yet):

Model comparisons:
elpd_diff se_diff
lModel 0.0 0.0
nlModel -190.7 17.9

I have been through the vignette and will try to go through it more thoroughly. Is there something in the model specification that is obviously wrong to you? All code for replicating the above is available here.

1 Like

Just a quick thing: the formula expects Angle to be in radians (between 0 and 2\pi). Will look into more details later. Best of luck!

1 Like

A quick update here. That seems like a big improvement to get closer to the data:
image

It still performs quite poorly when compared to the linear model using loo, though that may in part be explained by not taking “Part” into account. I can’t find an obvious way to address that but will work on that tomorrow. It could also be interesting to see if the model would perform better if we allow the angle and the centre to vary with the data.

Once again thanks for your help and it seems like we may indeed be very close (if not already there) to having solved the original problem.

It depends what you want to assume - either each particpant is allowed to have a slightly different shape of the ellipse, in which case you would have something like: a ~ 1 + 1|Part, b ~ 1 + 1|Part or you want the ellipse to be the same shape for everybody, but assume some people just have the ellipse larger or smaller, in which case you would have something like Distance ~ <original formula> + 1|Part.

Finally, since you are fitting strictly positive data, I would consider using lognormal family for both model variants.

Remember to post the final working code for posterity!

It seems like we are almost there, but the model does not converge every time in fact it more often than not fails to converge, stating high Rhats and low ESS.

I am quite certain that I have had the model converge at around a = 7 and b = 6 with both normal and lognormal priors:

Log normal priors

priors <- prior(lognormal(1, 1), nlpar = “a”) + prior(lognormal(1, 1), nlpar = “b”)

Normal Priors

priors <- c(set_prior(“normal(6.4,2)”, nlpar = “a”, lb = 0),
set_prior(“normal(6.4,2)”, nlpar = “b”, lb = 0))

One problem may be the assumption of A >= B which is not stated in the prior or model. I don’t know if that is possible?

The model is currently as follows:

brm(bf(Distance ~ sqrt(a ^ 2 * sin(AngleR)^2 + b^2 * cos(AngleR)^2),
a ~ 1+ 1|Part, b ~ 1+ 1|Part,
nl = TRUE), data = data, prior = priors, cores = 2, chains = 2, iter=6000, warmup = 2000,
control = list(adapt_delta = 0.999, max_treedepth = 15))

I have tried to find a seed that would reliably converge so I could share it but so far I have been out of luck.

Should I keep tweaking the priors or am I missing something obvious, here?

Will do.

Ahhh, the convergence problems. I thought that the parametrization might turn out to be problematic :-( I would be highly surprised if tweaking the priors helped, probably a reparametrization is needed. I am lazy to try to work out the math or run the model myself, but the thing to do is to inspect the pairwise posteriors (via pairs or using Shinystan) and look for suspicious patterns. (I describe some of the suspicious patterns at https://www.martinmodrak.cz/2018/05/14/identifying-non-identifiability/). The mcmc_parcoord could also be useful - example at https://mc-stan.org/bayesplot/articles/visual-mcmc-diagnostics.html. My best guess is that the model would work better parametrizing the ellipse via eccentricity, and some sort of “radius” but that might be wrong.

Also the lb = 0 might be also needed when using the lognormal prior (but maybe that’s added by brms as default, not sure.

AFAIK A >= B means the ellipse is wider than it is tall - why would you want to make such an assumption? You could however easily encode this by having parameter c with lower bound at 0 and substite a + c for all uses of b.

Hope that helps.