How to analyse this set of variables using causal modelling EDIT: edited the question slightly to keep the diagram intact but make it clear that we have an observed value and a hidden unmeasurable variable.
Here is the story. There are patients who have a certain disease D. We suspect that the disease might be aggraviated/caused by a certain condition C (our hypothesis: C causes D, $C \rightarrow D$). We cannot measure the condition C in the patients, however we can measure some variable $O$ which depends on $C$ ($C \rightarrow O$).
But herein lies the problem: patients who suffer from D receive (often) a medication $M$. This medication is correlated with the disease status, but it is not identical (there are patients without $D$ who also receive it, and vice versa, there are patients with $D$ who do not receive $M$). Therefore, we have the causal relationship $D \rightarrow M$.
Moreover, the medication $M$ might influence the observed variable $O$. At least that is what a reviewer might say: you are only observing an association between $D$ and $O$ because there is an association $D\rightarrow M$ and an independent association $M \rightarrow O$. This suggests the following, more complex model: we have the $C$ (which cannot be directly measured) that might cause $D$, and we have $O$ dependent on $C$, which might be influenced by $M$:
$$C \rightarrow D \rightarrow M \rightarrow O$$
$$C \rightarrow O$$
Now, we observe a clear association between $O$ and $D$. We also cannot rule out an association between $M$ and $O$. Furthermore, if I run a regression with $D$ as the response variable, and both $M$ and $O$ as predictors, I still get very similar results... but much, much weaker. That is, the p value is higher and the effect size smaller. This is not unexpected, because $M$ is strongly associated with $D$, and I am removing the portion of the variance in $D$ that can be explained by $M$.
To make the matter a bit more complicated, I have several of these $C$ variables, $C_1 ... C_n$ and corresponding observations $O_1 .. O_n$. In general, the p-values and effect sizes are correlated between the model which does not include $M$ as a predictor and the model that does include $M$ as predictor. Notably, the associations between $M$ and $O_1 ... O_n$ are all not significant when corrected for multiple testing (the associations between $D$ and $O_1 ... O_n$ are significant regardless).
My question is, how should I analyse and present the result of the analysis? What I did until now was to create these alternate models, with or without $M$ as predictor for all $O_1 ... O_n$, and also shown that the associations between $D$ and $O$ are present for each level of $M$. Again, this latter part is weaker firstly because the groups are smaller and secondly because the data is imbalanced (given that $D$ and $M$ are associated, there are relatively few patients who have not been medicated and few controls who have).
How do I convince myself (and others) that $C \rightarrow D$ and this effect is not due to $M \rightarrow O$? Would it be convincing enough to test O ~ M | D (i.e., check what is the association between $O$ and $M$ when accounting for D) and show that no results are significant? Is there a better way?
Another idea: how about testing O ~ M + D and comparing it to the model O ~ D? If there is a substantial effect $M \rightarrow O$, then these models should be significantly different? (I do not observe this).
Equivalently (I think), I could calculate $O'$ as the residuals of the model O ~ D and test whether O' ~ M.
What bothers me in these ideas is that I am searching for lack of significance, but absence of evidence is not evidence of absence, and hence I am looking for a more concrete solution.
P.S. To make the question more substantial, imagine that we build the following simulation in R:
set.seed(1234)
# hidden variable
C <- rnorm(100)

# disease depending on the hidden variable
D <- ifelse(C + rnorm(100) > 0, "y", "n")
D <- factor(D, levels=c("y", "n"))

# medication depends on disease
M <- ifelse(as.numeric(D) + .5 * rnorm(100) < 1.5, "y", "n")
M <- factor(M, levels=c("y", "n"))

# observed values
O <- C + .5 * rnorm(100)

# test whether D is associated with O
cor.test(O, as.numeric(D), method="s") # p < 1e-7

# test whether M is associated with O
cor.test(O, as.numeric(M), method="s") # p < 0.01

Now, the job is to show that there is an association $C\rightarrow D$ that is independent of $M$ using only $O$, $M$ and $D$ variables. We know from how the simulation is derived that $O$ does not depend directly on $M$, and the observed association is only due to the path $O \leftarrow C \rightarrow D \rightarrow M$.
 A: [EDIT] The following solution no longer applies, except for the argument that $M\not\to C.$ Original question has been edited. Not deleting because of valuable comments by Alexis.
I think time is a major factor, here. You are observing multiple variables at different times, and you have feedback. The regular tools of causality are challenging to work with in such a case. However, here's an option for you to try: model the same quantity at different times as essentially different variables. After all, you will have different values for them in your data. Here's the model I would propose in your slightly simpler setting of just one precondition. The parentheses here show the value of that variable at that time. The different between time $t$ and time $t+1$ can be whatever you want:

The edge $D(t)\to C(t+1)$ is debatable, but seems reasonable to me. You might also debate $C(t)\to C(t+1),$ but you did posit that that edge exists.
Now the basic rules of causality say that causes must precede effects. Moreover, $C(t+1)$ occurs later in time than $C(t).$ Since $C(t)\to D(t)\to M(t)\to C(t+1),$ it cannot be (by the laws of causality) that $M(t)\to C(t).$
As for the analysis, if you are interested in the causal effect $C(t)\to D(t),$ then you need only regress $D(t)\sim C(t).$ There is no backdoor path from $C(t)$ to $D(t),$ so long as you only include pre-medicated data on the $C$ variable.
