1

I am trying to fully understand and reproduce the paper Normative brain size variation and brain shape diversity in humans, in which the authors have used brain surface measurements from many different subjects with different brain sizes to compute a scaling factor for each vertex on the brain mesh that describes the increase/decrease in surface area at that anatomical location relative to a given increase in the whole brain surface.

In the Methods section, which is available in the Supplementary Materials, p4, the authors write:

Vertex-wise areal scaling with total cortical area was estimated using a semi-parametric spline regression model that simultaneously incorporated effects of age and sex on vertex area (41, 42). The scaling relationships presented in our results, were estimated using the following fixed-effects model:

Log10(vertex_area) ~ s(age, by=sex) + ß1[log10(total_area)], [1]

where s(age, by=sex) denotes a fixed-effect thin plate spline basis estimated for each sex. [...] The ß1 coefficient from model [1] provides > an estimate of how area at each cortical vertex scales with total cortical size.

I would like to implement that model (preferably in Python, Matlab or GNU R), but I am struggling with the part s(age, by=sex) denotes a fixed-effect thin plate spline basis estimated for each sex.

How could I compute such a fixed-effect thin plate spline basis, given observations of the position of the respective surface point for a number of subjects (and their age+sex)?

spirit
  • 113
  • 3

1 Answers1

1

The mgcv package in R has a gamm function which will provide what you need as it can incorporate so-called factor smooth interactions, which allow the smooth, nonlinear effect of one continuous variable such as age to be different across the levels of a categorical variable such as sex. See for example:

https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/smooth.construct.fs.smooth.spec.html

https://www.fromthebottomoftheheap.net/2017/10/10/difference-splines-i/

and especially

http://www.sfs.uni-tuebingen.de/~jvanrij/Tutorial/GAMM.html.

You'll probably need to specify your model along these lines:

gamm(log10(vertex_area) ~ s(age, by = sex, bs = "ts") +     
                          log10(total_area), 
                          data = YourData)

In fact, you could also use the bam() function of mgcv instead of the gamm() function, which is suitable for large data sets, as illustrated in the last of the 3 links I provided.

Isabella Ghement
  • 20,314
  • 2
  • 34
  • 58