Choosing covariate values for making predictions I have a model as follows
flexsurvreg(Surv(time, cens) ~ sex, data = "dist", data = data)

While plotting the model's prediction, it conforms well with Kaplan-Meier curves for both sexes.
When I run an adjusted model and do a prediction at mean age
flexsurvreg(Surv(time, cens) ~ sex + age, data = "dist", data = data)

The lines do not conform well at first months from follow up. The reason for this can be the different age distributions of sexes in the data. However, when using mean age - 3 years in predictions, lines from the second model and Kaplan-Meier conform very well.
Thus, should I publish:
a) the results on mean/median age and the predictions do not conform that well at early follow up with Kaplan-Meier curves
b) use slightly manipulated age for predictions, which seems to give more accurate predictions for both sexes? Or the manipulation leads to some biases?
Age ranges from 30-100 years in the data.
 A: Remember that your model is not fundamentally trying to predict the whole-population Kaplan-Meier curve. You are trying to estimate survival as functions of covariate values. Illustrate your model with sets of covariate values that make sense from your understanding of the subject matter. Don't fret that the corresponding curves don't match the Kaplan-Meier curves, as a Kaplan-Meier curve should not even be expected to agree with a parametric survival model evaluated at the mean value of a covariate.
One reason can be seen from the general form of a parameterized accelerated-life model:
$$\log T=-x'\beta + \sigma W $$
where $T$ is the time to event, $x'$ is the transpose of the vector of covariate values, $\beta$ is the vector of regression coefficients, $W$ is an error term from a specific distribution (e.g., standard normal for a log-normal model), and $\sigma$ is a scale factor. Thus a linear change in a covariate value leads to a corresponding change in the logarithm of an event time.
As a trivial example, if you have a single covariate with coefficient $\beta=1$ and two cases having covariate values of $x=1$ and $x=3$, at $W=0$ the average of their times to events $T$ is $(e^{-1}+e^{-3})/2 \approx 0.21$; at their mean covariate value of $x=2$, the corresponding time to event is $e^{-2} \approx 0.14$, a good deal shorter. To match the average between the two cases with $x=1$ and $x=3$, you would have to choose a covariate value of about $1.57$, lower than the mean of those covariate values. I suspect that something similar is going on with your data.
Furthermore, cases having covariate values associated with earlier event times will tend to be over-represented at the earlier times in Kaplan-Meier curves. Other things being equal, you will for example tend to have older individuals over-represented at early times in a typical overall-survival Kaplan-Meier curve. A major reason for semi-parametric Cox models or fully parametric survival models is to dissect apart how different covariates are associated with outcome, so that you can go beyond the simple overall summary of the study population provided by a Kaplan-Meier curve.
