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
Here across doctors, is an estimated covariance matrix, and the Gaussian errors have variance 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 . 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. The new model accounts for doctor-specific baselines and age slopes.
The command fits a Gaussian generalized additive mixed model
The population growth curve is a centered penalized regression spline with a cubic regression basis. The incubator term remains a random slope without a random intercept, because ~0+hours removes the intercept. A population linear component is allowed within the smooth's unpenalized part; curvature is penalized rather than being required.
A gamm fit returns a list with gam and lme components. Plot the population smooth and pointwise uncertainty bands with
plot(fly.model2$gam, se=TRUE, residuals=TRUE)
Inspect whether the fitted curve departs materially from an affine function, particularly in regions supported by observations. The effective degrees of freedom near one suggest that the penalization has selected an approximately linear centered smooth. Bands and partial residuals help distinguish convincing curvature from noisy deviations, although pointwise bands are not a simultaneous test of linearity. A formal comparison would need to account for both smoothing selection and the incubator dependence. The smooth plot tests the visual plausibility of a common straight growth curve while retaining incubator-specific slope variation.