I am trying to replicate the sleep deprivation study found in the R package LME4 (for fitting linear mixed effects models) and described in this paper. The dataset `sleepstudy`

contains

the average reaction time per day for subjects in a sleep deprivation study. On day 0 the subjects had their normal amount of sleep. Starting that night they were restricted to 3 hours of sleep per night. The observations represent the average reaction time on a series of tests given each day to each subject.

What I want to to is to fit a random effects model on intercept and slope for each subject. And I would like to use the matrix parametrization. That is, `X`

is the design matrix for the fixed effects given by the formula `Reaction ~ 1 + Days`

, and `Z`

is the design matrix for the random effects. However, I have some doubts on how to design such matrix `Z`

to get intercept’s and slope’s coefficients.

What I did is to simply create a design matrix from the formula `Reaction ~ 0 + Subject + Subject:Days`

. Is this the right way? On one hand the results I get seems reasonable (see notebook). On the other hand, in the paper describing LME4 (link) at page 9 they use a more complex method based on Khatri-Rao and Kronecker product. Am I missing something?