Can I estimate Variance of Gamma from Negative Binomial distribution distributed data, given NB is Gamma-Poisson compound

I believe the data I have follows Negative Binomial distribution (over-dispersed Poisson). We know Negative Binomial is a Gamma-Poisson compound distribution. The variance of this Gamma distribution is actually the variance I am interested in. Can I estimate the variance of the Gamma distribution of the compound from just my Negative Binomial data.

if $X|\lambda ~ \sim Poisson (\lambda)$

$\lambda \sim Gamma(shape = r, scale = p/(1 − p))$

then, the marginal distribution of X follows $NB(r,p)$,

$var(X) = \frac{pr}{(1-p)^2}$

$mean(X) = \frac{pr}{1-p}$

$var(\lambda) = r \cdot \frac{p}{1-p}$

Hence, $var(\lambda) = var(X) - mean(X) = \frac{pr}{(1-p)^2} - \frac{pr}{1-p} = \frac{rp}{1-p}$

Also, we can consider X and $\lambda$ are two random variables, then

$Var(X) = E(Var(X|\lambda))+ Var(E(X|\lambda))$

$Var(X) = E(\lambda)+ Var(\lambda)$

$Var(\lambda) = Var(X) - mean(\lambda)$

$E(X) = E(E(X|\lambda)) = E(\lambda)$

From the derivations, it seems like I can estimate the variance of Gamma distribution in the mixture by $Var(X)-mean(X)$.

Alternatively, if we parameterize Negative Binomial X using $\mu$ and size, the $var(X) = \mu + \mu^2/size$. So $var(\lambda) = \mu^2/size$.

The formulas you have given can be used to construct moment estimators. But you can also use maximum likelihood estimation in the marginal negative binomial distribution. The variance of the gamma distribution is a function of the negbin parameters, so just use the same function on the maximum likelihood estimates of $$r,p$$ to get the maximum likelihood estimator of the gamma variance (by invariance of maxlik estimators: Invariance property of maximum likelihood estimator?).