= Solution
There is a genuine syntax error in the printed command: `random=list(...)` must be a separate argument, following a comma rather than a plus. The intended call is
``
model2 <- gamm(npres ~ s(age, bs="cr") + s(wbc, bs="cr"), random=list(doctor=~age))
``
With the default Gaussian family, this <generalized additive mixed model> is
$$
\boxed{Y_{ij}=\alpha+s_1(a_{ij})+s_2(w_{ij})+b_{0j}+b_{1j}a_{ij}+\varepsilon_{ij}.}
$$
Here $(b_{0j},b_{1j})^T\overset{\rm iid}{\sim}N_2(0,D)$ across doctors, $D$ is an estimated <covariance matrix>, and the Gaussian errors have <variance> $\sigma^2$ and are independent of the <random effects>. The formula `~age` includes a <random intercept> as well as a <random slope>; it does not force these two effects to be uncorrelated.
The <random intercept> allows different prescribing baselines, and the <random slope> allows doctor-specific linear modifications to the population age curve. Conditional on these effects the observations are independent under the model, whereas two patients of the same doctor have additional <covariance> $(1,a_{ij})D(1,a_{kj})^T$. This accounts for <clustered data> and avoids treating doctor variation as independent patient-level noise. The population smooth and doctor-specific predictions are different targets. With only four doctors, estimates of the random-effect <covariance> can be imprecise; the Gaussian support issue for counts also remains. \b[The new model accounts for doctor-specific baselines and age slopes.]
Back to article page
