# Varying slopes model

In psychology, we increasingly encounter data that is nested. It is to the point now where any quantitative psychologist worth their salt must know how to analyze multilevel data. A common approach to multilevel modeling is the varying effects approach, where the relation between a predictor and an outcome variable is modeled both within clusters of data e.

And there is no better way to analyze this kind of data than with Bayesian statistics. Not only does Bayesian statistics give solutions that are directly interpretable in the language of probability, but Bayesian models can be infinitely more complex than Frequentist ones. This is crucial when dealing with multilevel models, which get complex quickly.

Stan is the lingua franca for programming Bayesian models. You code your model using the Stan language and then run the model using a data science language like R or Python. Stan is extremely powerful, but it is also intimidating even for an experienced programmer. I assume a basic grasp of Bayesian logic i. I also assume familiarity with R and the tidyverse packages in particular, ggplot2, dplyr, and purrr. This installation is more involved than typical R packages.

One more note before we dive in. Some time back I wrote up a demonstration using the brms package, which allows you to run Bayesian mixed models and more using familiar model syntax. I started with brms and am gradually building up competency in Stan.

Stan is the way to go if you want more control and a deeper understanding of your models, but maybe brms is a better place to start. We are creating repeated measures data where we have several days of observations for a group of participants each denoted by pid. Our outcome variable is y and our continuous predictor is x. This code will make more sense once we start working with the models.

We now have a data frame that is very much like something you would come across in psychology research. We can see a good amount of clustering at the participant level. This is not the way to analyze this data, but I use it as a simple demonstration of how to construct Stan code.

We then assign priors to the parameters. Priors encode our knowledge and uncertainty about the data before we run the model. True, the likelihood function i. But how do we know these are reasonable priors? We create a grid of evenly spaced values. We then want to see the likelihood of these x s in the context of three normal distributions with the same mean but different standard deviations.

## R-bloggers

We can see that the three distributions vary in terms of how spread out they are. We might say that the top distribution with mean 0 and sd 0. Flatter distributions allocate probability more evenly and are therefore more open to extreme values.This post builds on our recent introduction to multi-level modeling with tfprobability, the R wrapper to TensorFlow Probability.

In a previous postwe showed how to use tfprobability — the R interface to TensorFlow Probability — to build a multilevelor partial pooling model of tadpole survival in differently sized and thus, differing in inhabitant number tanks. A completely pooled model would have resulted in a global estimate of survival count, irrespective of tank, while an unpooled model would have learned to predict survival count for each tank separately.

The former approach does not take into account different circumstances; the latter does not make use of common information. Also, it clearly has no predictive use unless we want to make predictions for the very same entities we used to train the model. In contrast, a partially pooled model lets you make predictions for the familiar, as well as new entities: Just use the appropriate prior.

Assuming we are in fact interested in the same entities — why would we want to apply partial pooling? For the same reasons so much effort in machine learning goes into devising regularization mechanisms.

In the tadpole example, this means we expect generalization to work better for tanks with many inhabitants, compared to more solitary environments. For the latter ones, we better take a peek at survival rates from other tanks, to supplement the sparse, idiosyncratic information available.

Or using the technical term, in the latter case we hope for the model to shrink its estimates toward the overall mean more noticeably than in the former. This type of information sharing is already very useful, but it gets better. The tadpole model is a varying intercepts model, as McElreath calls it or random interceptsas it is sometimes — confusingly — called 1 — intercepts referring to the way we make predictions for entities here: tankswith no predictor variables present.

So if we can pool information about intercepts, why not pool information about slopes as well? This will allow us to, in addition, make use of relationships between variables learnt on different entities in the training set.

Unlike the tadpole case, this time we work with simulated data. For today, we stay with the simulated data for two reasons: First, the subject matter per se is non-trivial enough; and second, we want to keep careful track of what our model does, and whether its output is sufficiently close to the results McElreath obtained from Stan 2. So, the scenario is this. In terms of intercepts and slopes, we can picture the morning waits as intercepts, and the resultant afternoon waits as arising due to the slopes of the lines joining each morning and afternoon wait, respectively.

So is that all? Actually, no. In our scenario, intercepts and slopes are related. But first, we actually have to generate the data. Consider taking a look at the first section of that post for a quick reminder of the overall procedure. The same goes for the R packages tensorflow and tfprobability : Please install the respective development versions from github.

The first five distributions are priors. First, we have the prior for the correlation matrix.You can report issue about the content on this page here Want to share your content on R-bloggers? In psychology, we increasingly encounter data that is nested.

It is to the point now where any quantitative psychologist worth their salt must know how to analyze multilevel data. A common approach to multilevel modeling is the varying effects approach, where the relation between a predictor and an outcome variable is modeled both within clusters of data e. And there is no better way to analyze this kind of data than with Bayesian statistics. Not only does Bayesian statistics give solutions that are directly interpretable in the language of probability, but Bayesian models can be infinitely more complex than Frequentist ones.

This is crucial when dealing with multilevel models, which get complex quickly.

Stan is the lingua franca for programming Bayesian models. You code your model using the Stan language and then run the model using a data science language like R or Python. Stan is extremely powerful, but it is also intimidating even for an experienced programmer.

I assume a basic grasp of Bayesian logic i. I also assume familiarity with R and the tidyverse packages in particular, ggplot2, dplyr, and purrr. This installation is more involved than typical R packages. One more note before we dive in. Some time back I wrote up a demonstration using the brms package, which allows you to run Bayesian mixed models and more using familiar model syntax.

I started with brms and am gradually building up competency in Stan. Stan is the way to go if you want more control and a deeper understanding of your models, but maybe brms is a better place to start. We are creating repeated measures data where we have several days of observations for a group of participants each denoted by pid. Our outcome variable is y and our continuous predictor is x. This code will make more sense once we start working with the models. We now have a data frame that is very much like something you would come across in psychology research.

We can see a good amount of clustering at the participant level. This is not the way to analyze this data, but I use it as a simple demonstration of how to construct Stan code. We then assign priors to the parameters. Priors encode our knowledge and uncertainty about the data before we run the model.

True, the likelihood function i. But how do we know these are reasonable priors? We create a grid of evenly spaced values. We then want to see the likelihood of these x s in the context of three normal distributions with the same mean but different standard deviations.

We can see that the three distributions vary in terms of how spread out they are. We might say that the top distribution with mean 0 and sd 0. Flatter distributions allocate probability more evenly and are therefore more open to extreme values. We need to look at what kind of data is compatible with the priors. Remember, all you need to create a line is an intercept and a slope!

The trick this time is to generate intercepts and slopes from the different normal distributions. The more restrictive set of priors on the top constrains the lines to have intercepts and slopes close to zero, and you can see that the lines vary little from one another. This lets you take advantage of autocompletion and syntax error detection. The data block is where we define our observed variables.We have seen how random intercept models allow us to include explanatory variables and we saw that, just like with the variance components model, in the random intercept model, each group has a line, and we saw that the group lines all have the same slope as the overall regression line.

And remember that was true for the variance components model as well, because in that case all the lines were flat, they just had slope 0. So, for the random intercept model, in every group, the effect of the explanatory variable on the response is the same and that's actually one of the assumptions of the random intercept model.

Well here is a possible situation we could have, we've got some data points here, and we can imagine that these are exam results for pupils within schools, so along the x axis, we have previous exam score, and along the y axis we have exam score at age 16 and we want to fit a random intercept model to this data.

So here's our random intercept model and we fit it to our data and in order to better be able to examine how well that model fits the data, we're just going to highlight four of the groups and look at those. So if we look at the red group, here, you can see that for this group, the points are following a line with a steeper slope than the group line that we've drawn in, and again for the dark blue group, the points seem to be following a line with a steeper slope than the group line that we've drawn in.

On the other hand for the light blue group, the points seem to be following a line with a shallower slope than the group line that we've drawn in. And for the green group as well, the points seem to be following a line with a shallower slope than the group line that we've drawn in.

So for this data, for some groups, the explanatory variable has a large effect on the response and for others it has a small effect. So clearly the random intercepts model, with its parallel group lines, is not doing a very good job of fitting the data.

## Random slope models

So that's all very well in theory, but you might wonder, does this actually happen in practice? Well, in fact, in exactly the example we've been considering, pupils within schools with response being exam score, and explanatory variable being previous exam score, some investigators have found that their data behaves like this.

So that for some schools pretest has a large effect on the response and for others the effect is smaller. But on the other hand, other investigators with exactly the same situation, pupils within schools; response: exam score; explanatory variable: previous exam score, have found that for their data, the random intercepts model is a perfectly good fit: it doesn't appear that the relation between the explanatory variable and the response is different for different schools.

And also it's important to bear in mind that for some datasets, there's only enough power to fit a random intercepts model in any case. So what this tells us is that sometimes the random intercept model does fit well to the data and we don't want to look any further, that's perfectly adequate, but in other cases the random intercept model doesn't fit well and we need something else. So for the data we've just been looking at, what we really want is a model that looks like this.

So this is a random slopes model and we can colour in the other groups as well.

### Hierarchical partial pooling, continued: Varying slopes models with TensorFlow Probability

You can see that this is actually fitting the data better than the random intercept model did because you can see for example for the red group, the group line now does seem to have the same slope as the line that the points are following, and for the light blue group as well, the group line seems to have the same slope as the line that the points are following, and the same for all the other groups.

And we can extend those group lines as well, the model doesn't specify how long those lines are, they extend infinitely, but, as you can see, drawing it like that does look a bit complicated and confusing, it's kind of hard to see what's going on there. So we mostly will not extend the lines when we draw it on future slides. Well, unlike a random intercept model, a random slope model allows each group line to have a different slope and that means that the random slope model allows the explanatory variable to have a different effect for each group.

It allows the relationship between the explanatory variable and the response to be different for each group. So how do we actually achieve that in terms of the model equation? Well what we do, is to add a random term to the coefficient of x 1 so that it can be different for each group. For the random slope model we have added in this u 1 x 1This u 1 is different for every group, so that means that this coefficient is different for every group.

So that means that the relationship between x 1 and y is different for every group. So, what do our u 1 and our u 0 look like in our random slopes model?

Well, our u 0 looks just the same as for the random intercept model, it's still the difference between the intercept for the overall line and the intercept for the group line, so just this difference here. So, to be able to be able to better see what u 1 is we've actually drawn faintly here a line parallel to the overall line.

It's not really part of the model, this, we've just drawn it on so it's easier to see what u 1 is. So the slope for the group line is this, so u 1 is the difference between the slope for the group line and the slope for the overall line. If we look at the light blue group now, again u 0 is the difference between the intercept for the overall line and the intercept for the group line, this difference here, and now for this group, the group line has a shallower slope than the overall line.In many cross-national comparative studies, mixed effects models are being used in which a number of slopes are fixed and the slopes of one or two variables of interested are allowed to vary across countries.

The aim is often then to explain the varying slopes by referring to some country-level characteristic. My question is whether it is possible that the estimation of these varying slopes the interesting ones is affected by the fact that the slopes of other uninteresting variables are fixed even though they may actually vary over countries and if so, how it may be affected? Do you think the decision not to include many varying slopes is predominantly methodologically informed?

And do you think Bayesian analyses can provide a better solution for the kind of situations where many slopes should be allowed to vary? I do not think the reviewer is correct, but I would not know how to react to this. In response to your first question: essentially what you have is an omitted variable bias so, yes, if you have varying coefficients that you treat as being fixed, then this can cause problems.

Second, when many slopes vary, it can be necessary to do more regularization. The flat priors we use in my book with Jennifer might not be enough.

I think we need more research on informative priors for hyperparameters such as, in this case, group-level covariance matrices. Most of the models I use have varying slopes for all variables as it is the default choice when using Gaussian processes.

Can you give an example with a minimum working example with code, I mean? The main usage for Gaussian Processes is to describe distributions over functions.

Suppose for example you have brightness of a signal across a row of pixels. Then each pixel could be at a brightness that is essentially independent of the other pixels. You put a prior on L that makes it about the same size as the pixel array… and this tells the Gaussian process machinery that you expect things to change slowly from one end of the array to the other.

Such a long length scale smooth covariance function puts high probability on functions that look like low-order polynomials, say 2nd to 5th order polynomials over the data points. There is no hard cutoff that prevents something like a 6th or 8th order polynomial from coming out of samples of this gaussian process, but the majority of these samples will be well approximated by lower order polynomials. More specifically, gaussian processes are meaningful to me when there is a continuous manifold over which some independent variable can vary… like for example x and y positions of pixels, or latitude and longitude, or time.

But it does matter. If you have wild estimates, the SEs of your fixed effects are going to be wider. Especially in low-power situations pretty much the norm at least in psycholinguisticsyou are going to miss effects at an alarming rate. Of course, having more data lets you see smaller effects with more clarity, but the importance of that will vary a lot from research project to research project. Sometimes there is something going on that can be seen even with small data sets!

Your suggestion makes a lot of sense in general, but would not work in areas like psychology and psycholinguistics.Covariance matrices allow us to capture parameter correlations in multivariate hierarchical models; sampling these using Hamiltonian Monte Carlo in Tensorflow Probability can be tricky and confusing; this post is about some of the math involved and how to get this right.

Hierarchical models allow us to account for variations between different groups in our data. Varying intercepts models allow us to fit different models to different tanks, while pooling together information between tanks. Varying intercepts are already very powerful models. However, in many most?

R - Multilevel Model Example

Each model has two parameters a slope and an interceptand we allow these to vary. We can also allow them to covary. For example, if higher slopes usually go with lower intercepts, we want to know that, and use that to improve our estimation of both. Luckily for us, inLewandowski, Kurowicka, and Joe published a method for generating random correlation matrices, aptly referred to as the LKJ distribution.

So far so good - sampling correlation matrices seems straightforward. The problem starts when we want to use a Hamiltonian Monte Carlo and we usually want to use Hamiltonian Monte Carlo to sample from some larger model that contains an LKJ distribution.

HMC allows us to generate samples from arbitrary joint distributions, not only from distributions for which we have explicit sampling methods. We add a small helper function to avoid rewriting all the kernels everytime; see the previous post for explanations about the different function calls here.

Without getting into the gory mathematical details 2a bijector can be thought of as a differentiable one-to-one mapping between the unconstrained space in which the HMC trajectories live, and the constrained manifold. HMC produces samples in the unconstrained space, and the appropriate bijector spits out a valid correlation matrix.

For us, this bijector is tfb. Cholesky factors come up in many different places in statistics, machine learning, metric learning, computational linear algebra, etc. We generate samples from a bivariate guassian with zero mean, unit variance and no correlation:.

We now define a correlation matrix between two variables with correlation We compute its Cholesky factor, multiply it with the original uncorrelated data, and voila:.

So the bijector solves the constrained-unconstrained problem, and HMC can run smoothly. The HMC sampler works with the log probability function of the model. But LKJ is a distribution over correlation matrices, not Cholesky factors of correlation matrices, which is the output of our bijector! So what can we do? To do so, we need two more bijectors: tfb.

Chain which, surprisingly, chains composes the two bijectors:. This looks good. But this is cumbersome, and not very readable. So we end up sampling cholesky factors, tranforming them back to correlation matrices just to compute their cholesky factors again…. Since tfp-nightly It also helps in verifying the model is correctly specified and that the MCMC sampler does what you think it does, which is good.

The average morning waiting time is the intercept, and the difference between afternoon and morning average waiting times is the slope.Below is a hypothetical dataset with 15 students nested in 3 schools. Their passing an exam is correlated with their family income and the school they attend to which is the nesting variable. If I want to write a regression analysis with c onstant intercept and slope for each school:. I would:. If I want to write a regression analysis with varying intercept and slope for each school :.

View solution in original post. No, you wouldn't do this unless you had a very good reason. You want all terms in one model, which gives you better estimate of the overall variability than you get if you do it with a BY statement. I don't see this as answering the original question. This discussion points to the root of my original problem: I think I'm having a bit of a hard time to understand the difference between having a "varying slope" vs.

Ksharpthanks for your response as well! If you are interested in the entire population of schools and you have randomly selected these schools and you want to know the variability of the intercept across schools, or the variability of slopes across schools in other words, a standard deviation of the intercepts or a standard deviation of the slopesthen you have a RANDOM effect.

Need further help from the community? Please sign in and ask a new question. Learn about the latest articles, videos, and code samples to keep your SAS skills fresh! Sign In. Turn on suggestions. Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for. Search instead for. Did you mean:. All forum topics Previous Next. Varying intercept and slope in regression analyses.

Posted PM views. Thanks a lot in advance!

Accepted Solutions. Re: Varying intercept and slope in regression analyses. Reeza wrote: 2 - BY statement. Thanks a lot for your response Paige! Yes, that's what I meant -- Paige Miller. Ksharp wrote: that wil lead you to Generalize Mixed Model. Thanks a lot in advance to both of you Sign in. Recommended by SAS. For personalized recommendations, sign in with your SAS profile. Discussion stats.