# What are the sampling distributions underlying a product of gaussian variables (moderated effect)?

I'm trying for a while now to understand the differential power of main effects (two random variables) and moderated effects (the product of these two random variables). At this point, I'm mostly interested in Gaussian. I'm working with a model like : $$$$y=\beta_0 +\beta_1 x+\beta_2 z+\beta_3 xz+\epsilon_y$$$$ where $$$$\begin{bmatrix} x \\ z \end{bmatrix} \sim \mathcal{N} \left( \begin{bmatrix}\mu_X \\ \mu_Z \end{bmatrix} ,\,\begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix} \right)$$$$ I'm specifically looking for the sampling distribution of $$\beta_3$$ compared to regression coefficients (like $$\beta_1$$). Acknowledging that $$t$$-distributions should be used to test their power, I have found that there are many differences. Because $$t$$ is a normal ($$\mathcal{N}$$) divided by $$\sqrt{\chi^2/n}$$, I investigated in this direction, and found out these differences.

First, the normal distribution of $$\beta_3$$ as slightly different standard deviation than $$\beta_1$$, i.e., $$\sqrt{\frac{\sigma^2}{(n-8)(1+\rho^2)}}$$ (or something very similar) instead of $$\sqrt{ \frac{ \sigma^2}{n}}$$.

Second, the expected $$\chi^2$$ seems to be a lognormal instead. It has heavier tails and is more skewed. I have not found how to derive the parameters of the log normal yet, but I've come across some approximations.

Finally, the $$t$$-distribution of $$\beta_3$$ does not behave as expected. It has a slightly underestimated noncentral parameters ($$\delta$$) and, more importantly, its variance is not the expected $$v/(v-2)(1+\delta^2)-\mu_1^2$$, $$v$$ being the degree of freedom and $$\mu_1$$ is its first moment, except when $$\delta=0$$. Using the previous points and this answer, I have tried to fit the function:

$$f(z)=\frac{1}{2\pi\sigma_x\sigma_y}\int_0^\infty \exp\{ -[zy-\mu_x]^2/2\sigma_x^2-[\log(y)-\mu_y]^2/2\sigma_y^2\}~\text dy$$

where $$Z$$ is the normal, $$Y$$ the lognormal, and dividing $$\mu_y$$ and $$\sigma_y$$ by 2, because I am interested in the square root of the lognormal (see this question). However, this pdf does not yield "perfect" results (see last figure).

My questions are :

1. Is there a way to analytically derived the parameters of the lognormal based on the informations in the models?
2. What would be the function to fit the "$$t$$-distribution" of $$\beta_3$$ or how to correct the graphic presentation, or is it normal to be that way, or is the distribution actually a skewed generalized $$t$$ distribution?

This problems is very important to me and any help would be appreciated. Thank you,

As a supplement, this is an example of code on what I'm trying to achieve.

library(foreach)
library(doParallel)
library(fitdistrplus)
cores = detectCores()
cl = makeCluster(cores[1]-1)
registerDoParallel(cl)

# The function in https://stats.stackexchange.com/a/507175/102655
# Z = normal
# X = lognormal
f = function(z,x,mz,mx,sz,sx) {
exp(- (z * x - mz)^2/(2 * sz ^ 2)-(log(x) - mx) ^ 2 / (2 * sx ^ 2))
}
norm_lognorm = function(z,mz,mx,sz,sx){
1/(2*pi*sz*sx) * integrate(f = f, lower=0, upper=Inf,
z = z, mz = mz, mx = mx, sz = sz, sx = sx,
rel.tol=.Machine$$double.eps^.5)$$value
}

# Parameters
reps = 100000 # number of replication
n = 50        # sample size
b = .50       # the parameters of interest (beta_3)
rho = .5      # correlation between main effects

txz = foreach(l = 1:reps, .combine = cbind) %dopar% {
x = rnorm(n)
z = rho * x + sqrt(1 - rho ^ 2) * rnorm(n)
e = rnorm(n)
y1 = b * x * z + sqrt(1 - b ^ 2*(1 + rho ^ 2)) * e
txz = summary(lm(y1 ~ x * z))$$coefficient[4,1:3] w = rnorm(n) y2 = w * b + sqrt(1 - b ^ 2) * e tw = summary(lm(y2 ~ x + z + w))$$coefficient[4,1:3]
txz = c(txz, tw)
}
# txz contains estimate of xz, standard error of xz, t-value or xz, ...
# estimate of w, standard error of w, t-value of w.
# xz is the product of two random gaussian; w is a single gaussian variable to compare

Beta = matrix(c(0,0,b))
Sw   = matrix(c(1,rho,0,rho,1,0,0,0,1),3,3)
Sxz  = matrix(c(1,rho,0,rho,1,0,0,0,1 + rho^2),3,3)

# R2 for both models
R2w  = c(t(Beta)%*%Sw%*%Beta)
R2xz = c(t(Beta)%*%Sxz%*%Beta)

# estimation of the parameters of the lognormal
est = summary(fitdist(txz[2,]^2, distr = "lnorm", method = "mle"))$estimate # Would like to know them a priori, even though I have found approximation sx = est[2] # mx = est[1] # # parameters for the normal mz = b sz = sqrt(1 - R2xz) / (sqrt((n - 8)*(1 + rho^2))) # graph par(mfcol = c(3,1), mai = c(.25, .25, .25, .25),oma = c(4, 4, 1, 1)) # blue is for xz and red if for w hxz = hist(txz[1,], plot=FALSE, breaks=50) # hist of estimates of xz, blue hw = hist(txz[4,], plot=FALSE, breaks=50) # hist of estimates of w, ref m = c(min(hxz$$breaks, hw$$breaks), max(hxz$$breaks, hw$$breaks)) # min and max of x-axis plot(hxz, freq=FALSE, main=NULL, ylim = (c(0,max(hxz$$density,hw$$density))), xlim = m, col=rgb(0,0,1,.1)) plot(hw, freq=FALSE, add=TRUE, col = rgb(1,0,0,.1)) xabs = seq(m[1], m[2], by=sum(abs(m))/10000) lines(xabs, dnorm(xabs, mean = b, sd = sqrt(1-R2w) / (sqrt(n))), col = "red",lty=1) lines(xabs, dnorm(xabs, mean = mz, sd = sz), col="blue", lty=1) legend("topright", legend = ("\u03b2"), box.lty = 0, text.font = 3, bg = "transparent") hxz = hist(txz[2,]^2, plot=FALSE, breaks=75) # hist of variances of est. of xz, blue hw = hist(txz[5,]^2, plot=FALSE, breaks=50) # hist of variances of est. of w, red m = c(min(hxz$$breaks,hw$$breaks), max(hxz$$breaks,hw$$breaks)) plot(hxz, freq=FALSE, main=NULL, ylim = (c(0,max(hxz$$density,hw$$density))), xlim = m, col = rgb(0,0,1,.1)) plot(hw, freq = FALSE, add = TRUE, col = rgb(1,0,0,.1)) abx = seq(m[1], m[2], by = sum(abs(m))/10000) y = dchisq((abx)*(n-4)*(n-1)/(1-R2w), df=(n-4)) lines(abx, y*(n-4)*(n-1)/(1-R2w), lty=1, col="red") y = dlnorm(abx, meanlog = mx, sdlog = sx) lines(abx, y, lty=1, col="blue") legend("topright", legend = ("var(\u03b2)"), box.lty = 0, text.font = 3, bg = "transparent") hxz= hist(txz[3,], plot=FALSE, breaks=50) # hist of t values of est. of xz, blue hw = hist(txz[6,], plot=FALSE, breaks=75) # hist of t values of est. of w, red m = c(min(hxz$$breaks,hw$$breaks), max(hxz$$breaks,hw$$breaks)) plot(hxz, freq=FALSE, main=NULL, ylim = c(0, max(hxz$$density,hw$$density)), xlim = m, col = rgb(0,0,1,.1)) plot(hw, freq = FALSE, add = TRUE, col = rgb(1,0,0,.1)) abx = seq(m[1], m[2], by = sum(abs(m))/10000) yy = dt(abx, ncp = b*sqrt(n-4)/sqrt(1-R2w), df = n-2) lines(abx, yy, col = "red") yy = apply(as.matrix(abx), MARGIN=1, FUN = norm_lognorm, mz = mz, mx = mx/2, sz = sz, sx = sx/2) lines(abx, yy, col="blue") legend("topright", legend = ("t"), box.lty = 0, text.font = 3, bg = "transparent")  The blue is related to $$\beta_3$$ and the red to $$\beta_1$$. The simulation use $$n=50, \rho=.5, \beta_1=\beta_3=.5$$ (both $$\beta$$ are used in different data sets). I am confident for the results in red histograms and red lines, but not necessarily for the blue ones. Thanks for your helps! • That code is ugly. You could make it easier to read by explaining what$b$is, by replacing the repeated stuff with functions of$ha$and$hb\$, and by adding some spaces (eg as recommended at adv-r.had.co.nz/Style.html). Jul 10 at 12:55
• Thank for the input, I have updated my code. I hope it is clearer and cleaner!
– POC
Jul 11 at 1:02
• That's better, but the definition of norm_lognorm is still hard to read, many individual lines still can only be read with left-right scrolling, and the code for the graphics still has a lot of repetition -- which can respectively be improved by using dnorm, more line breaks, and more functions. With clean code, there might be a clean answer; as it stands now, it's likely that some simple error is obscured by the mess, just like the typo in "stackexchxznge". Jul 11 at 3:27
• Hope it is better!
– POC
Jul 11 at 3:36
• Hi @Ben! Yes, I think I am looking at the sampling distribution of the OLS estimator.
– POC
Jul 13 at 1:29