Model-averaged estimates/intervals/distributions

[strengejacke]

cf #768

[strengejacke]

library(MATA)
library(parameters)

set.seed(0)
n = 20                              # 'n' is assumed to be even
x1 = c(rep(0,n/2), rep(1,n/2))      # two groups: x1=0, and x1=1
x2 = rnorm(n, mean=10, sd=3)
y = rnorm(n, mean = 3*x1 + 0.1*x2)  # data generation

x1 = factor(x1)
m1 = glm(y ~ x1)                    # using 'glm' provides AIC values.
m2 = glm(y ~ x1 + x2)               # using 'lm' doesn't.
aic = c(m1$aic, m2$aic)
delta.aic = aic - min(aic)
model.weights = exp(-0.5*delta.aic) / sum(exp(-0.5*delta.aic))
residual.dfs = c(m1$df.residual, m2$df.residual)

p1 = predict(m1, se = TRUE, newdata = list(x1 = factor(0), x2 = 10.64), type = "link")
p2 = predict(m2, se = TRUE, newdata = list(x1 = factor(0), x2 = 10.64), type = "link")
theta.hats = c(p1$fit, p2$fit)
se.theta.hats = c(p1$se.fit, p2$se.fit)

#  95% MATA-Wald confidence interval for theta:
mata.wald(theta.hats, se.theta.hats, model.weights, mata.t = TRUE, residual.dfs = residual.dfs)
#> [1] 0.6196654 1.6516695

parameters:::averaged_parameters(m1, m2)
#> [1] 0.619622 1.651641

^{Created on 2022-09-05 with reprex v2.0.2}

[strengejacke]

@bwiernik implementation seems to work, but I'm not sure about the meaning of this function. It returns one estimate (in the above example currently only CIs)? What does this estimate represent?

[codecov-commenter]

Codecov Report

Merging #771 (d44c8c7) into main (4e62652) will decrease coverage by 0.13%.
The diff coverage is 0.00%.

@@            Coverage Diff             @@
##             main     #771      +/-   ##
==========================================
- Coverage   52.99%   52.85%   -0.14%     
==========================================
  Files         184      185       +1     
  Lines       12706    12738      +32     
==========================================
  Hits         6733     6733              
- Misses       5973     6005      +32     

Impacted Files	Coverage Δ
R/averaged_parameters.R	0.00% <0.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[bwiernik]

No, that's not right. It should be output with the same structure as the input vectors, averaged across models (eg, if 32 predictions from 2 models, then output should be 32 predictions, weighted across the two models). The MATA CIs for the 32 predictions would then be the values for each prediction that satisfy the critical value when averaged across models

[github-actions]

[air] _{reported by reviewdog 🐶}

Suggested change

	f = .tailarea, interval = c(-1e+10, 1e+10), theta_hats = theta_hats,
	se_theta_hats = se_theta_hats, model_weights = model_weights, alpha = alpha,
	residual_dfs = residual_dfs, tol = 1e-10
	f = .tailarea,
	interval = c(-1e+10, 1e+10),
	theta_hats = theta_hats,
	se_theta_hats = se_theta_hats,
	model_weights = model_weights,
	alpha = alpha,
	residual_dfs = residual_dfs,
	tol = 1e-10

[github-actions]

[air] _{reported by reviewdog 🐶}

Suggested change

	f = .tailarea, interval = c(-1e+10, 1e+10), theta_hats = theta_hats,
	se_theta_hats = se_theta_hats, model_weights = model_weights, alpha = 1 - alpha,
	residual_dfs = residual_dfs, tol = 1e-10
	f = .tailarea,
	interval = c(-1e+10, 1e+10),
	theta_hats = theta_hats,
	se_theta_hats = se_theta_hats,
	model_weights = model_weights,
	alpha = 1 - alpha,
	residual_dfs = residual_dfs,
	tol = 1e-10

[github-actions]

[air] _{reported by reviewdog 🐶}

Suggested change

	.tailarea <- function(theta, theta_hats, se_theta_hats, model_weights, alpha, residual_dfs) {
	.tailarea <- function(
	theta,
	theta_hats,
	se_theta_hats,
	model_weights,
	alpha,
	residual_dfs
	) {