Your models 1-3 ignore one or both main effects (except when your variables are coded as factors; see https://stackoverflow.com/questions/40729701/how-to-use-formula-in-r-to-exclude-main-effect-but-retain-interaction). You rarely want to ignore main effects. Since only model 4 contains all main effects and the interaction, this is most likely the most appropriate one. This is usually written as
For a discussion of interaction effects without main effects, look for example here: Including the interaction but not the main effects in a model.
Edit: I actually just learned that
lme4 would not even let you estimate a random slope if you have just two time points, because you are trying to estimate just as many random effects as you have observations (or even fewer if there is missing data). As Ben Bolker points out, although with
nlme the estimation appears to work, the model will not be able to actually distinguish between random slope variance and residual variation.
With health interventions whose outcome develops over time, the standard is to use an ANCOVA adjusted for baseline, which will turn out to be quite similar to the mixed model with random intercept. So if you have measured at baseline (i.e. at the time of randomization, or when the treatment was assigned), the usual way to estimate the intervention effect would be:
lm(variable ~ baseline + time*group)
Or if you are really interested in the variance of the random intercept, you could use
lme(variable ~ baseline + time*group, random= ~ time|subject)
but that will give (almost) identical results as the other model, and the variance of the random intercept will be (almost) the same as the variance of the baseline values.
Be also aware that the interpretation of the results can be somewhat challenging in the presence of interaction effects. You might want to read further into the topic, and examine the fitted model further with pairwise comparison (e.g. with the package
emmeans) to get a better understanding of the group- and time-specific estimates.