Isolating Effects in a Causal Model I am working on trying to recover correct group means in a causal inference model (I'm pretty new to causal inference) that I'm running on a simulated dataset. I believe that the issue I'm running into is the method and order in which I am controlling for different effects in the dataset.
Broadly speaking, here's the model that I'm interested in fitting.

We have data divided into 4 groups. For each observed value of obs, we get four pieces of information: which individual individual 1 was, which individual individual 2 was, and the value of the covariate effect. Each individual in the collections of individual 1 and individual 2 has some unknown, unobservable effect mean; the values of  effect are observed for each observation, as is which group it came from. Additionally, note that the means of the effect values are associated with an individual 1 value, but have some random fluctuation.
The means of the unobserved values of individual 1 and individual 2 vary by group, as does the mean of the values of effect. We're interested in recovering the group means of individual 1 after controling for individual 2 and effect.
Here is the data generating process in R:
set.seed(20)
group <- c(rep(1, 100), rep(2, 100), rep(3, 100), rep(4, 100))
individual1 <- 1:400
individual2 <- 1:400
rand_effect1 <- c(rnorm(100, mean = 0, sd = 1),
                  rnorm(100, mean = 0.5, sd = 1),
                  rnorm(100, mean = 1, sd = 1),
                  rnorm(100, mean = 1.5, sd = 1))
rand_effect2 <- c(rnorm(100, mean = 0, sd = 1),
                  rnorm(100, mean = -1, sd = 1),
                  rnorm(100, mean = -2, sd = 1),
                  rnorm(100, mean = -3, sd = 1))
fix_effect <- c(rnorm(100, mean = 0, sd = 1),
                rnorm(100, mean = 0.25, sd = 1),
                rnorm(100, mean = 0.5, sd = 1),
                rnorm(100, mean = 0.75, sd = 1))
ind1 <- sample(individual1, size = 5000, replace = TRUE)
re1 <- rand_effect1[ind]
fe <- fix_effect[ind] + rnorm(5000, mean = 0, sd = 0.5)
grp <- group[ind]
f <- (floor((ind1-1) / 100) * 100) + 1
c <- (floor((ind1-1) / 100) + 1) * 100
ind2 <- c()
for(i in 1:5000) {
  v <- floor(i / 100)
  ind2 <- c(ind2, 
            sample(f[i]:c[i], 1))
}
re2 <- rand_effect2[ind2]
obs <- re1 + re2 + fe + rnorm(5000, mean = 0, sd = 4)
training_data <- data.frame(
  group = as.character(grp),
  individual1 = as.character(ind1),
  individual2 = as.character(ind2),
  effect = fe,
  obs = obs
)

I could certainly throw all of this into a random effects model with individual 1 and individual 2 handled appropriately. However, it would be ideal to include group means in the model; this model should be able to handle individuals with small sample sizes from different groups appropriately. Knowing nothing about an individual 1 value except for the fact that they came from group 3 should be able to tell us, eventually, that their individual random effect is in the ballpark of 1 (as is clear in the data generating process). However, it's clearly not as easy as throwing in a group factor into the model, since there are different group effects on all three inputs.
Any suggestions about how to appropriately control for each component in this model would be very appreciated!


However, it's clearly not as easy as throwing in a group factor into the model, since there are different group effects on all three inputs.

I see only one way how group has an effect and that is in the different means of the normal distributions.
rand_effect1 <- c(rnorm(100, mean = 0, sd = 1),
                  rnorm(100, mean = 0.5, sd = 1),
                  rnorm(100, mean = 1, sd = 1),
                  rnorm(100, mean = 1.5, sd = 1))
rand_effect2 <- c(rnorm(100, mean = 0, sd = 1),
                  rnorm(100, mean = -1, sd = 1),
                  rnorm(100, mean = -2, sd = 1),
                  rnorm(100, mean = -3, sd = 1))
fix_effect <- c(rnorm(100, mean = 0, sd = 1),
                rnorm(100, mean = 0.25, sd = 1),
                rnorm(100, mean = 0.5, sd = 1),
                rnorm(100, mean = 0.75, sd = 1))

You can just as well consider normal distributions with zero mean and a separate fixed effect for group.
rand_effect1 <- c(rnorm(100, mean = 0, sd = 1),
                  rnorm(100, mean = 0, sd = 1) + 0.5,
                  rnorm(100, mean = 0, sd = 1) + 1,
                  rnorm(100, mean = 0, sd = 1) + 1.5)
rand_effect2 <- c(rnorm(100, mean = 0, sd = 1),
                  rnorm(100, mean = 0, sd = 1) -1 ,
                  rnorm(100, mean = 0, sd = 1) - 2,
                  rnorm(100, mean = 0, sd = 1) - 3)
fix_effect <- c(rnorm(100, mean = 0, sd = 1),
                rnorm(100, mean = 0, sd = 1) + 0.25,
                rnorm(100, mean = 0, sd = 1) + 0.5,
                rnorm(100, mean = 0, sd = 1) + 0.75)

Or is this just a simplification of your true problem and do you have more complicated ways how group has an effect?

In addition, because the three different effects of group all result in effectively a single effect, you can not obtain the 12 individual means. You can only discover the sum of three means and obtain 4 parameters. Or if you mix individuals of different groups then you should be able to figure out 8 parameters.
