Is it possible to use a latent variable calculated from a fitted regression as the dependent variable in another regression?

Hello,

I have a logged linear regression that predicts total muscle “U” from a measurement of length “L” and adipose fat “A”.

\log(U) \sim \alpha_{lu} + \beta_{lu1} \ast \log(L) + \beta_{lu2} \ast \log(A).

Once fitted, this equation can be used to predict the mass of structural muscle (Ku) which is associated with structure (bones etc) as:

\log(Ku) \leftarrow \alpha_{lu} + \beta_{lu1} \ast \log(L)

From which I can estimate the amount of storage muscle (Ou) which is associated with adipose mass (energy stores)

Ou \leftarrow U - Ku

The proportion of storage that is muscle (PrOou) can then be found as:

PrO_{Ou} \leftarrow \frac{Ou}{Ou+A}

I would like to then use this calculated, latent variable to fit a beta regression that relates a measure of body condition (BC) (overall fatness) to the proportion of storage that is muscle (PrOou)

PrO_{Ou} \sim beta(Au, Bu)
logit(\mu) \leftarrow \alpha_{pr} + \beta_{pr} \ast BC
Au = (\mu)*phi_{U}
Bu = (1 - \mu))*phi_{U}

I have tried to make this work multiple ways but I don’t seem to be able to use the PrOou value calculated in any of the blocks to be the dependent variable in the beta regression.

The model below runs but although PrOou is calculated in the transformed parameter block it does not appear that the beta regression uses it, instead it appears to be fit to nothing or missing data as the estimated parameters for the beta regression are completely uniformed.

Any advice or discussion on how I might best fit these “sequential” regressions would be much appreciated! Even if I can figure out how to use the latent variable as a dependent in the beta regression I wonder if the likelihood could be sorted out or if I might be better to just take a Generalized likelihood uncertainty estimation approach.


//multiple linear regression code for Logged total muscle regression

data {
  int<lower=0> nbears;   // number of dissected bears
  int<lower=0> KlU;   // number of total muscle linear regression predictors
  
  matrix[nbears, KlU] lUPreds;   // total muscle linear regression predictor matrix where each row is an observation and each column a predictor
  vector[nbears] lU;      // Observed logged total muscle
  
  vector[nbears] lL; // Observed logged straight line body length predictor for lKnU
  vector[nbears] lKnU;   // Observed logged total structure non-muscle
  
  vector[nbears] U; //Observed mass total muscle
  vector[nbears] A; //Observed mass of adipose
  
  vector[nbears] SMI; //Scaled mass index calculated from observed total mass and length
}


parameters {
  real alphalU;           // total muscle linear regression intercept
  vector[KlU] betalU;      // total muscle linear regression coefficients
  real<lower=0> sigmalU;  // total muscle linear regression error scale
  
  real alphalKnU;  // total structure non-muscle intercept
  real betalKnU; //total structure non-muscle slope with logged length
  real<lower=0> sigmalKnU;  // total structure non-muscle linear regression error scale
  
  real alphaPrU; //Proportion of storage that is muscle intercept
  real<upper=0> betaPrU; //Proportion of storage that is muscle slope with SMI
  real<lower=0, upper=300> phiU;                // Proportion of storage that is muscle dispersion parameter
  
}

transformed parameters{
  
  vector[nbears] Ku; //estimate for the proportion of storage that is muscle
  vector[nbears] PrOu;
  
  Ku= exp(alphalU + betalU[1]*lL); //Estimate the mass of structural muscle 
  PrOu = (U-Ku) ./ (U-Ku+A); //Estimate the proportion of storage that is muscle

}

model {
    // model calculations
  //vector[nbears] Xbeta;                  // linear predictor
  //vector[nbears] mu;                     // transformed linear predictor
  vector[nbears] Au;             // shape parameter for beta distn
  vector[nbears] Bu;             // shape parameter for beta distn
   
  // priors
  phiU ~ uniform(0, 300);          // put upper on phi if using this
  
  // likelihood
    
    lU ~ normal(lUPreds* betalU + alphalU, sigmalU);  // logged total muscle likelihood
  
    lKnU ~ normal(betalKnU * lL + alphalKnU, sigmalKnU);  // logged total structure non-muscle likelihood

    for (i in 1:nbears){
      
      Au[i] = (inv_logit(alphaPrU + SMI[i] * betaPrU))*phiU;
      Bu[i] = (1-inv_logit(alphaPrU + SMI[i] * betaPrU))*phiU;
  
      PrOu[i]~beta(Au[i], Bu[i]);
      }
      
}

// track variables of interest
generated quantities {
  
      vector[nbears] muKu; 
      vector[nbears] muPrOu; 
  
      muKu= exp(alphalU + betalU[1]*lL); //Estimate the mass of structural muscle  
      muPrOu = (U-muKu) ./ (U-muKu+A); //Estimate the proportion of storage that is muscle
}

1 Like

So, this bit:

sets the value of the PrOu vector as a transform involving both data variables U and A and the parameter Ku.

Then this bit:

achieves a prior on PrOu.

So if this is what you want, there’s already an issue: Since you PrOu is a transformed parameter and it’s transform is non-linear with respect to the parameter Ku, to place a prior on PrOu you’ll have to add a jacobian adjustment.

Next, if you want PrOu to be constrained by data in some way, you need it to be on the right-side of a sampling statement that has data on the left-hand side.

1 Like

Thank you for your reply. I didn’t realize it was setting a prior which is great to know.

My problem is that the data to constrain PrOu can only be calculated with the fitted variables from:

lU ~ normal(lUPreds* betalU + alphalU, sigmalU);

So rather than fitting this equation, then estimating PrOu and imputing it as data into a new model to fit as a beta regression, I’m trying to find a way to estimate PrOu within the script and fit it to a beta regression at the same time.

The goal is to have a fitted beta regression that can predict PrOu from SMI along with the fitted log regressions that predict lU and lKnU

UPDATE:

I’ve now tried a new approach to do this using a “dummy” place holder in the data{} block where I feed in PrOu=1 for all individuals. I then calculated PrOuT as a transformed parameter to modify this data placeholder.

I multiply the dummy data PrOu, a vector of ones, by the calculated value for PrOuT on the left hand side of the beta regression to attempt to fit it.

This way seems to “work” in terms of recognizing (PrOu[i]*PrOuT[i])~beta(Au, Bu) as an equation to be fit rather than a prior BUT it now gives me unrealistic values for the parameters in the fitted lU regression.

Perhaps I can still use this trick but need to add a jacobian adjustment ?

//multiple linear regression code for Logged total muscle regression
data {
  int<lower=0> nbears;   // number of dissected bears
  int<lower=0> KlU;   // number of total muscle linear regression predictors
  
  matrix[nbears, KlU] lUPreds;   // total muscle linear regression predictor matrix where each row is an observation and each column a predictor
  vector[nbears] lU;      // Observed logged total muscle
  
  vector[nbears] lL; // Observed logged straight line body length predictor for lKnU
  vector[nbears] lKnU;   // Observed logged total structure non-muscle
  
  vector[nbears] U; //Observed mass total muscle
  vector[nbears] A; //Observed mass of adipose
  
  vector[nbears] SMI; //Scaled mass index calculated from observed total mass and length
  
  vector[nbears] PrOu; //DUMMY VARIABLE, vector of ones to multiply the transformed parameter estimate by
}


parameters {
  real alphalU;           // total muscle linear regression intercept
  vector[KlU] betalU;      // total muscle linear regression coefficients
  real<lower=0> sigmalU;  // total muscle linear regression error scale
  
  real alphalKnU;  // total structure non-muscle intercept
  real betalKnU; //total structure non-muscle slope with logged length
  real<lower=0> sigmalKnU;  // total structure non-muscle linear regression error scale
  
  real alphaPrU; //Proportion of storage that is muscle intercept
  real<upper=0> betaPrU; //Proportion of storage that is muscle slope with SMI
  real<lower=0, upper=300> phiU;                // Proportion of storage that is muscle dispersion parameter
  
}

transformed parameters{
  
  vector[nbears] Ku; //estimate for the proportion of storage that is muscle
  vector[nbears] PrOuT; //parameter to modify the DUMMY VARIABLE by
  
  Ku= exp(alphalU + betalU[1]*lL); //Estimate the mass of structural muscle using the fitted parameter from lU regression 
  PrOuT = (U-Ku) ./ (U-Ku+A); //Estimate of the proportion of storage that is muscle

}

model {
    // model calculations
  vector[nbears] Au;             // 1st shape parameter to relate beta distribution to linear regression 
  vector[nbears] Bu;             // 2nd shape parameter to relate beta distribution to linear regression 
  
   
  // priors
  phiU ~ uniform(0, 300);          // put upper on phi if using this
  
  
// likelihood
    
    lU ~ normal(lUPreds* betalU + alphalU, sigmalU);  // logged total muscle likelihood
  
    lKnU ~ normal(betalKnU * lL + alphalKnU, sigmalKnU);  // logged total structure non-muscle likelihood
      

      Au = (inv_logit(alphaPrU + SMI * betaPrU))*phiU;
      Bu = (1 - inv_logit(alphaPrU + SMI * betaPrU))*phiU;
      
      for (i in 1:nbears){
        (PrOuT[i]*PrOu[i]) ~ beta(Au, Bu);
      }
    
}

// track variables of interest
generated quantities {
  
      vector[nbears] muKu; 
      vector[nbears] muPrOu; 
  
      muKu= exp(alphalU + betalU[1]*lL); //Estimate the mass of structural muscle  
      muPrOu = (U-muKu) ./ (U-muKu+A); //Estimate the proportion of storage that is muscle
}

I think there’s still some confusion here. What do you mean by “beta regression”?

Essentially you use a link function (logit in my case here) to fit a linear regression to the beta distribution.

\mu = inv_logit (\alpha_{PrU} + SMI * \beta_{PrU})

which is transformed into the two beta shape parameters using the precision \phi

Au= \mu*\phi
Bu = (1-\mu)*\phi

https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.13234

" Beta regression is a technique that has been proposed for modelling of data for which the observations are limited to the open interval (0, 1) (Ferrari & Cribari-Neto, 2004; Smithson & Verkuilen, [2006]"… .“Beta regression consists of the same three components as generalized linear models (GLMs) (Bolker et al., [2009] McCullagh & Nelder, [1989], and those familiar with GLM will recognize the most important aspects of beta regression (the distinction between the two arises from the non-orthogonality of the model parameters, see below). Here, we briefly review these three elements: the random component (the beta distribution and its implied mean–variance relationship), the systematic component (the linear predictor) and the link function (specifying the link between the random and systematic component)…”

Here’s some examples in STAN
https://www.rpubs.com/kaz_yos/stan_beta1

Thank you!

It’s more complicated than this, both because because the transform is a many-to-one function of parameters, and because even if a Jacobian were available, applying the Jacobian would interfere with the priors on the terms in the first regression (which are currently improper, but might get contorted into something quite weird if the only prior is the beta regression plus a Jacobian adjustment).

@staned what this means is that you are not going to be able to get meaningful inference from this model without thinking very carefully about what priors you actually want. If you can come up with a self-consistent prior, then we can figure out how to write down the model.

I also worry that the math here seems to be a bit weird. Substituting your second equation into your first and rewriting for clarity, we have \log(U) = \log(K) + \beta_2\log(A) + \epsilon , where \epsilon \sim N(0, \sigma). Then you additively decompose U into a structural part and an adipose part:

It sorta looks like you are envisioning a decomposition where there’s a structural part that depends on log(K) and then an adipose part that corresponds to \beta_2\log(A) + \epsilon. But these two parts are additive on the log scale, which means that they’re multiplicative on the identify scale. So your estimate of O has a large dependency on K and not just on \beta_2\log(A) + \epsilon. Is that really what you’re going for here?

1 Like

Ah, I need to find a way to remember that such transforms add a wrinkle. I also didn’t fully understand the beta regression stuff, hence my initial thinking that the PrOu[i]~beta(Au[i], Bu[i]); bit was intended as a prior of some sort. Thanks @jsocolar for lending your greater expertise!

1 Like

hah, yeah right

I’m looking to parameterize a body composition model. So total mass is the sum of structural mass and storage mass. I have dissection data and everything but muscle can be readily assigned to storage or structure, but muscle is found in both storage and structure.

Essentially I’m trying to say that the observed total muscle is the sum of the unobserved, latent variables Ku (structural muscle) and Ou (storage muscle).

The first regression tries to describe total muscle based on predictors of structural mass (length L) and storage mass (adipose mass A) so that I can estimate the mass of structural muscle by setting A=0 in the fitted regression. I do this on the logged scale as they should be additive whereas on the normal scale they are multiplicative because a longer individual is more likely to have more adipose stores, we see a change in how these variables interact across values.

Your understanding of the decomposition is correct but I believe the adipose part would correspond to setting L=0 so would be Ou= exp(\alpha_{lU} + \beta_{lU2}*\log(A)) ?

The problem is that in regularly collected data I will not have the measurement of Adipose mass, only total mass and length so I cannot predict Ou using the logged muscle equation.

So given the estimate of structural muscle (Ku) I estimate what storage muscle must be given total muscle, leading to the estimate of the proportion of storage that is muscle (PrOu).

Ku= exp(\alpha_{lU} + \beta_{lU1}* \log(L)); //Estimate the mass of structural muscle using the fitted parameter from lU regression
PrOuT = \frac{(U-Ku)}{(U-Ku+A)}; //Estimate of the proportion of storage that is muscle

BUT I also won’t have total muscle mass is regular data so using that’s why I also need a way to predict our best guess of PrOu using a regression fit to a metric of body condition, SMI, which can be calculated using body length and total mass, measures that don’t require sacrificing an animal.

To summarize, I need to be able to predict Ku and KnU from length (L) only, and the estimate of PrOu from SMI (which is body condition based on L and total mass) to be able to have a body composition model that can function with non-destructive data.

Given that all these parts of the body sum to total mass there will be large dependencies in all observed variables, which is why I didn’t want to simply fit the first regression independent of the beta regression.

Thank you both very much for giving me the opportunity to explain further and for your willingness to think about this problem!

I am reading up on self-consistent priors at your suggestion but have not heard of this before. I do have a very good idea of what the parameter estimates in all regressions here should be based on other mammalian systems but I’m not sure if that’s helpful to develop the self consistent prior.

1 Like

What I understand is that there’s some total amount of muscle T that decomposes into a structural part S and a fat-storing part F. That is, T = S + F.

This implies to me that you don’t think adipose mass should be a strong predictor of structural mass. But because your regression is on the log scale, you are effectively multiplying your structural mass terms by your storage mass terms to get the total mass. This seems not to correspond with your conceptual model. If you want these terms to be additive, then you must add them on the identity scale. At a first brush, this might look something like:
T = e^{\alpha} + e^{\beta_1\log(L)} + e^{\beta_2\log(A)} + e^\epsilon
But then you point out that:

I don’t know if you mean that in longer individuals F should be higher given A, or in longer individuals A should be higher in general. In the former case, you can modify the third term in T = e^{\alpha} + e^{\beta_1\log(L)} + e^{\beta_2\log(A)} + e^\epsilon above to be e^{\beta_2\log(A)*\log(L)} or e^{\beta_2\log(A) + \beta_3\log(A)*\log(L)} as you see fit. In the latter case, you don’t need to change anything, as the effect is already captured in your measurement of A.

But there’s one final significant conceptual issue here. When we write T = e^{\alpha} + e^{\beta_1\log(L)} + e^{\beta_2\log(A)} + e^\epsilon, the second term corresponds to structural muscle, and the third term (or something like it; see above) corresponds to fat-storing muscle. But how do we partition the intercept term and the error term?

The intercept cannot be wholly both structural and storage. Likewise the error term cannot be wholly both structural and storage. You’ll need to decide how to partition the intercept to one, the other, or some combination, and likewise the error term.

Okay so basically using a GLM log link rather than logging both sides. I had thought that because variance would be expected to increases with an increase in T that I should use the Log-Log framework.

I think adding the interaction term would definitely work, it was something I was exploring before going log-log.

I thought that by setting A=0 (or 0.1 to avoid log(0)), the fitted equation would then only include the effect of structure and thus give the estimate for muscle that is associated with structure. In reality if length is 0 and adipose is 0 then we would have no muscle at all, would this be a case in which I should set the intercept to 0?