What is Hierarchical Linear Modeling?

This is going to be a very brief conceptual discussion on hierarchical linear modelling (HLM), also known as multi-level modelling in social sciences. People sometimes call this mixed-effects modelling or mixed modelling, but be careful that not all mixed-effects models have a hierarchy. 

Note: Mixed-effects models by definition are models which have both random and fixed effects (see below).

(1)  When do we use it? 

In social sciences, the data types are usually hierarchical in nature, i.e. they have a certain hierarchy that is based on grouping or nesting. We use HLM precisely in cases where observations are nested or clustered in some ways at a different level. It has many applications in cross-sectional and longitudinal research (repeated-measures design). Research studies create data hierarchies through the way the data are sampled. For example:

• Students (S) are nested within a particular Class, which again is nested in a School.
• Time points are nested within participants in a longitudinal study.

Fig-1: Pictorial representation of (a) two-level hierarchy and (b) three-level hierarchy (Source: SAGE Ency. 2018)

Suppose we wanted to examine A-level exam scores this year to predict the university entrance test in the following semester. So entrance exam is the outcome variable [DV] and the A-level score is the predictor [IV]. One way is to randomly sample students from the whole population such that each individual is guaranteed to be independent. However, what if the students recruited only belong to a particular high school? In HLM, the students are treated as the first level or the base, and the school as the second level (2-level hierarchy). In other words, we performed sampling primarily from four different local schools only. To advance further, we have to introduce 'Class' in between both levels (now, 3-level hierarchy), since each school surely has different classes. 

In HLM, independence can no longer be freely assumed. Why? E.g. think of how students of Teacher John of School ABC have a different profile, intelligence, general ability than those attending the class of Teacher Mary from School XYZ. In such a nested condition, assumptions fulfilled for the OLS will be violated: observations are no longer independent of each other, and residual errors are correlated within clusters. The error variance will be different for different clusters. As a result, the standard errors of the regression coefficients will generally be underestimated and the significance level will be wrong, leading to misinterpretation. One way is to bring the students' data of Level 1 up to the next level through some averaging, but this is a poor practice. HM treats levels as something to take into account. It is basically an extension of the ordinary least square (OLS) of the linear regression.   

(2)  How to construct the model? 

For a given A-level exam score x_i in a simple linear regression, we estimate the intercept b₀ and slope b₁. The term e_i is the residual. The equation states that for a unit increase in the A-level exam score, the score on the entrance test y_i increases by the value of b₁. But wait! The outcome dataset (entrance test scores y_i) is clustered by the school. If this is not taken into account, then the data from all students is treated as unique observations.  

Suppose we adopt a 2-level HLM comprising Students and Schools, i.e. the students are grouped or nested within a school. We introduce a group-level variation by changing the suffices, such that y_ij is the score on the entrance test for student i in school j, and u_j is the group-level residual for that j-th school. This u_j is a new random variable, referred to as the Level-2 or group residuals, and is assumed to follow a normal distribution, ~ N(0, σ²). The new HLM equation has the following terms:

• We see b₀ as the grand mean of the university entrance test y.  
• Level 1 contains the basic form of a linear regression, y_ij = b_0j + b₁x_ij + e_ij  ; in turns, Level 2 consists of b_0j = b₀ + u_j.
• The mean of y for the school group j is b₀ + u_j, where u_j is the "school effect", i.e. the difference between the mean of the school group j and the grand mean. With this, each school can have different mean scores on the entrance test. 
• The individual-level residual e_ij is the difference between the value of y for a student i and the individuals group mean b₀ + u_j. This reflects differences in students’ individual test scores from their respective school (or group) mean.

What can this improved model tell us? Some schools will have means that are higher than the grand mean, suggesting that their students perform better on the A-level exam on average, and some schools will have lower cluster mean values. 

(3)  Defining fixed and random terms 

The core of hierarchical models is the assignment of fixed effects and random effects, two terms that already appeared in another blog post on group-level fMRI analysis. You have to define each term of a mixed-effects model to be either a fixed or random effect.

• A fixed effect is a parameter that is fixed across all groups and does not vary in the model. It is the term of interest in experimental manipulation. Estimating a fixed-effect of a term is like estimating a regression slope b. Differences in the slope If all terms in the model above are fixed, it behaves as the usual linear regression.
• In contrast, a random effect allows each group to have a different estimate. In HLM, random effects represent a higher level variable under which data points are grouped. This implies that random effects must be categorical (but cannot be continuous!). For example, the residual errors (u_j and e_ij) are considered random effects. Estimating a random effect is like looking for an effect in our data to come from a large group of normally distributed datasets. 

The most basic model has only an intercept without any predictors x_ij, which is also called the intercept-only model. Another name for it is the Unconditional Model or Null Model. Usually, the model has only participants e_ij as a random error.

The next basic model is the one shown in the equation above. It allows for different groups or schools j to have a varying intercept or mean b_0j = b₀ + u_j. Thus, the model is also known as the random intercept model. 

(4)  The more complete model 

A random intercept model shown above assumes that the relationship between the university test and the A-level exam score (the predictor) is the same for each group. This assumption can be relaxed by allowing for different slopes for the predictor in each group, making it a random slope model as shown in the right panel of Fig-2 below. Each slope will be estimated separately for each group. In this new equation, there is a new term u_1j x_ij , where:

• The intercept for school group j is now b₀ + u_0j. 
• The slope for school j is now b₁ + u_ij , where b₁ is the average slope across groups.
• Subscript “0” differentiates the random effect for the intercept u_0j from the random effect for the slope u_1j. Both random effects are assumed to be normally distributed. 

  Fig-2: Hierarchical models [Source: Univ Bristol]

The variance of the intercept and slope are assumed to be correlated; the covariance between the intercept and slope is estimated as part of the random slopes model. The random slopes model is also commonly known as the random coefficient model or a growth curve model when using repeated measures or longitudinal analysis. 

(5)  Final notes 

How is the implementation in practice? Apart from defining fixed/random terms, we have to know which one nested under which variable.

• When we deal with a dataset we have to begin with the simplest model, which is the intercept-only model.
• Then we add additional terms, the simplest being the random intercept model, i.e. each group has its own group mean. 
• If we have more than one predictor or IV, define the fixed/random effect. A mixture of fixed and random effect slopes is possible. 
• The most complete model is the random slope model, where intercepts and slopes are treated as random effects. Often in a longitudinal dataset, we can add the quadratic term "Time" (Level-1).

So how do we know which model is best? We can use the usual performance metrics such as Bayesian Information Criterion (BIC) or Akaike Information Criterion (AIC). The two most common libraries for HLM analysis in R are lme4 and nlme. In the commands, you have to indicate which ones are the fixed and random effects, and which variable is nested within what. The outputs of the HLM in R is quite similar to what we expect from OLS in linear regression. 

The random slope model often mimics reality, but this requires a large sample size. Parameter estimation in HLM is achieved using either the Maximum Likelihood (ML) or Restricted Maximum Likelihood (ReML). ReML works better if the sample size is small.

Some good references:

[1]  https://rpsychologist.com/r-guide-longitudinal-lme-lmer

[2]  https://www.lcampanelli.org/mixed-effects-modeling-lme4/

[3]  https://m-clark.github.io/mixed-models-with-R/introduction.html

Introduction to Artificial Neural Network (1)

Everyone is talking about AI and machine learning recently. But the concepts of AI goes back to more than 40 years ago when scientists took steps to model the brain and cognition. In its original concept, there lies a perceptron, an artificial neuron that can form a neural network to learn and make decision (after F. Rosenblatt). Just as a neuron has an all-or-none firing characteristic, a neural network is able to provide a binary response or yes/no. While in the past it required advanced C++ and expensive hardware to run modeling, today knowing Python and having a high-speed computer are sufficient.

1. Principles of a perceptron
A perceptron takes several binary inputs and produces a single binary output. How does it compute the output? Each input is associated with a certain weight w_i and the total value of w_i x_i will be compared against a certain threshold to decide the output. Look at the figure below.

We can take a simplified real life example for how a perceptron works. Suppose there is a jazz festival and you need to decide whether you want to go ('1') or not ('0'). There seems to be 3 important factors that influence your decision, each has its own influence value:
      - Is the weather good (yes/no)? Influence scale = 10.
      - Is your girlfriend going with you (yes/no)? Influence scale = 3.
      - Is the place near a metro station (yes/no)? Influence scale = 4.
Note that the influence scales here represent the weights you personally assign to these factors. For example, you really hate being caught in the rain so much; or you don't really mind if you go alone without your girlfriend and make new friends. You should also infer a threshold to make your decision. That's it, a perceptron has just been used to model how you make a decision. What happens when we combine multiple perceptrons? This structure is called a neural network, as shown in the figure below. In practice, such network is used to take into account more factors and make a better decision.

Let us describe the perceptron more formally. First, the summation of weights and inputs can be rewritten as a dot product of w.x = Σ w_i x_i . Second, the threshold is in fact a bias, b = −threshold. This bias can be thought as a measure of how easy it is to get the perceptron to give an output '1'. Or to put it in more biological terms, it is a measure of how easy it is to get the perceptron to fire. The formula to get the output thus becomes, y = w . x + b. Thirdly, the binary output of '0' or '1' beyond a certain threshold mimics a type of function called a step function.

2. Simple neural network
A neural network is basically a multi-layer perceptron or MLP as shown in the figure above. The leftmost layer is called the input layer and the neurons inside it are called input neurons. The rightmost layer is called an output layer with just a single output neuron. In between, there exists 2 layers of what is called hidden layers, which take inputs from the input layer and project outputs to the output layer. A neural network in such structure is also popularly known as a feedforward neural network since the outputs from one layer are projected to the the next layer to the right and so on, with no feedback allowed. There are a few other neural network architecture, which will be said at a later time.

The power of a neural network lies in the capability to be trained. Neural networks can be 'trained' to behave in a certain way, or to produce a certain output given some patterns of inputs. Engineers can call this an iterative fine-tuning process. We use neural networks to help us decide something. For example, you want to decide whether a blurred image is a picture of a handwritten digit (number "0" to "9"). Training the neural network here means we let the network to readjust the weights and biases within the network according to the given inputs, in this case, different images.

How can we readjust the weights and biases? We let a small change in a weight or bias to cause only a small change in output Δy (like fine-tuning). In this way, we would get our network to behave more in the manner we want. For example, suppose the network was mistakenly classifying an image as an "8" when it should be a "9". We could figure out how to make a small change in the weights and biases so the network gets a little closer to classifying the image as a "9". This process is then repeated over and over again to make the output more and more accurate to decide "9". The network is said to be learning.

But there is one problem. With the current setup, a small change in an input x_i will yield quite a big jump in the binary output, it can completely be flipped (yes/no, '1' or '0'). This is primarily due to the all-or-none nature of the perceptron output. Scientists were not satisfied with such characteristic, so they defined a new type of perceptron where the output follows a sigmoid rather than a step function characteristic. Thus, a neuron with a sigmoid function don't just produce output '1' or '0'. It turns out that with this characteristic, Δy is a linear function of the changes Δw and Δb. This linearity makes it easy to choose small changes in the weights and biases to achieve any desired small change in the output. Formally, output characteristic of an artificial neuron or perceptron can be defined by the activation function, where a sigmoid is one of the examples.

As mentioned, training a neural network means allowing the weights w and biases b to readjust themselves such that a desired output is obtained. We can track how well the training progresses through a metric called a cost function, sometimes also a loss or objective function as a function of weights and biases. Remember that the goal of a neural network is to help us make decision or prediction. Traditionally, the cost function defines the difference between the predicted and the actual output of the network. Our end goal is to find w and b such that this difference is minimized.

3. Minimizing cost function
Suppose we have our multi-layer perceptron and this is so-called our model. Mathematically, the method to "train" the model of a neural network is called backpropagation algorithm. This term is coned after Rumelhart et al who proposed an efficient numerical solution of the training problem. Training or learning in terms of backpropagation here means to minimize discrepancy or error (or something bad) by adjusting weights and biases. 

To portray the discrepancy or error, we need a cost function, which defines how 'good' our model is at the moment as learning progresses. Now we wish to obtain a set of parameters (weights and biases) such that the discrepancy between the predicted and actual output of the model (neural network) as defined by the cost function is minimized. To minimize this so-called cost function, we can use an optimization algorithm called gradient descent, by iteratively moving in the direction of steepest descent as defined by the negative of the gradient (or slope).

Suppose we have a cost function F that depends on imaginary parameter x (x₁ and x₂). Let us imagine a valley defined in the dimension of x₁ and x₂, such as the one shown below, and we want to roll a ball down this valley. Making the ball roll down at different dimension of x₁ and x₂ is akin to saying that the change in F is negative, ΔF < 0. We can build a relationship such that: ΔF ≈ ∇F ⋅ Δx ; where ∇F is a the gradient vector of F, which carries partial differentiation operators. At the moment, it is enough to think a gradient vector simply as something that relates changes in x to changes in F, just as we would expect something called a gradient to do. Suppose we make the change in x to be Δx = −η∇F, where η > 0. Then mathematically, ΔF ≈ −η ∇F⋅∇F = −η∥∇F∥² , which means that ΔF is guaranteed to be negative and F will forever decrease not increase. The ball is for sure rolling down the valley!

The whole idea of iteration is as follows. From an arbitrary ball position of in dimension x, we first compute the change in Δx. so that to find a new position of x → x' = x − η∇F. Note that the arrow denotes an update rule, the variable takes up a new value. This update rule can be thought as defining the gradient descent algorithm. It gives us a way of repeatedly changing the ball position in order to find a minimum value of the function. If we keep doing this over and over again, we will keep decreasing F, until theoretically we reach a global minimum. Once this global minimum has been found, it is said that the training algorithm has converged. Note that η is called learning rate where it is usually kept small, to control learning and to prevent the update behavior to be chaotic.

To summarize, the way the gradient descent algorithm works is to repeatedly compute the gradient vector ∇F, and then to move in the opposite direction step by step, so as to "fall down" the valley. Updates of weight and bias parameters occur in an efficient way until the error is minimized.

Some notes:

1) In more complex networks for a multiclass prediction, softmax is used as the activation function instead of a sigmoid.

2) The following website is informative to understand backpropagation.