Translating R lme comand to mathematical equation I would appreciate if someone could help me in translating the following R command into a mathematical equation:
lme(score ~ factor(timeslot), random=~1|subjectid, data=a)

For each subject there are observations at 6 different times (timeslots) during a day.
 A: $$
y = X\beta + Z\gamma + \epsilon
$$
$y$ is an $N$ length vector of all observations across all subjects
$X$ is an $N \times 6$ design matrix. The first column is all $1$s for the intercept term which is implicit in the call to lme - i.e. your call is equivalent to
lme(score  ~ 1+factor(timeslot),random=~1|subjectid,data=a)
The remaining $5$ columns of $X$ map the observations to the timeslots at which they were observed. Note that there are $5$ rather than $6$ of these columns as the effects of each timeslot are taken relative to the first one.
$\beta$ is a $6$ length vector of fixed-effects where the first element is the intercept term and the remaining $5$ elements correspond to the relative effect of the last $5$ timeslots compared to the first.
$Z$ is an $N \times n$ design matrix mapping the observations to the subject on which they were measured. $n$ is the total number of subjects in your data set (which you have not specified).
$\gamma$ is an $n$ length vector of random intercept terms, one for each subject.
$\epsilon$ is an $N$ length vector of error terms.
