Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Modify model priors when using include.always #41

Merged
merged 2 commits into from
Nov 5, 2019
Merged

Modify model priors when using include.always #41

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a2a3982
adjusted model priors for nuisance predictors
vandenman Apr 17, 2019
abb73a6
added unit test
vandenman Apr 17, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/bas.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,8 @@ double trunc_beta_binomial(int modeldim, int p, double *hyper);
double trunc_poisson(int modeldim, int p, double *hyper);
double trunc_power_prior(int modeldim, int p, double *hyper);
double Bernoulli(int *model, int p, double *hyper);
double compute_prior_probs(int *model, int modeldim, int p, SEXP modelprior);
int no_prior_inclusion_is_1(int p, double *probs);
double compute_prior_probs(int *model, int modeldim, int p, SEXP modelprior, int noInclusionIs1);
void compute_margprobs_old(Bit **models, SEXP Rmodelprobs, double *margprobs, int k, int p);
void compute_modelprobs(SEXP modelprobs, SEXP logmarg, SEXP priorprobs, int k);
void set_bits(char *bits, int subset, int *pattern, int *position, int n);
Expand Down
3 changes: 2 additions & 1 deletion src/glm_deterministic.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@ SEXP glm_deterministic(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
int *model = (int *) R_alloc(p, sizeof(int));
memset(model, 0, p*sizeof(int));

int noInclusionIs1 = no_prior_inclusion_is_1(p, probs);
k = topk(models, probs, k, vars, n, p);

/* now fit all top k models */
Expand All @@ -82,7 +83,7 @@ SEXP glm_deterministic(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
SEXP glm_fit = PROTECT(glm_FitModel(X, Y, Rmodel_m, Roffset, Rweights,
glmfamily, Rcontrol, Rlaplace,
betapriorfamily));
double prior_m = compute_prior_probs(model,pmodel,p, modelprior);
double prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);
logmargy = REAL(getListElement(getListElement(glm_fit, "lpy"),"lpY"))[0];
shrinkage_m = REAL(getListElement(getListElement(glm_fit, "lpy"),
"shrinkage"))[0];
Expand Down
5 changes: 3 additions & 2 deletions src/glm_mcmc.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ SEXP glm_mcmc(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
n = sortvars(vars, probs, p);
for (i =n; i <p; i++) REAL(MCMCprobs)[vars[i].index] = probs[vars[i].index];
for (i =0; i <n; i++) REAL(MCMCprobs)[vars[i].index] = 0.0;
int noInclusionIs1 = no_prior_inclusion_is_1(p, probs);

// fill in the sure things
int *model = ivecalloc(p);
Expand Down Expand Up @@ -80,7 +81,7 @@ SEXP glm_mcmc(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
SEXP glm_fit = PROTECT(glm_FitModel(X, Y, Rmodel_m, Roffset, Rweights,
glmfamily, Rcontrol, Rlaplace,
betapriorfamily));
prior_m = compute_prior_probs(model,pmodel,p, modelprior);
prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);

logmargy = REAL(getListElement(getListElement(glm_fit, "lpy"),"lpY"))[0];
shrinkage_m = REAL(getListElement(getListElement(glm_fit, "lpy"),
Expand Down Expand Up @@ -135,7 +136,7 @@ SEXP glm_mcmc(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
glm_fit = PROTECT(glm_FitModel(X, Y, Rmodel_m, Roffset, Rweights,
glmfamily, Rcontrol, Rlaplace,
betapriorfamily));
prior_m = compute_prior_probs(model,pmodel,p, modelprior);
prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);

logmargy = REAL(getListElement(getListElement(glm_fit, "lpy"),"lpY"))[0];
shrinkage_m = REAL(getListElement(getListElement(glm_fit, "lpy"),
Expand Down
7 changes: 4 additions & 3 deletions src/glm_mcmcbas.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ SEXP glm_mcmcbas(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
n = sortvars(vars, probs, p);
for (i =n; i <p; i++) REAL(MCMCprobs)[vars[i].index] = probs[vars[i].index];
for (i =0; i <n; i++) REAL(MCMCprobs)[vars[i].index] = 0.0;
int noInclusionIs1 = no_prior_inclusion_is_1(p, probs);

// fill in the sure things
int *model = ivecalloc(p);
Expand Down Expand Up @@ -82,7 +83,7 @@ SEXP glm_mcmcbas(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
SEXP glm_fit = PROTECT(glm_FitModel(X, Y, Rmodel_m, Roffset, Rweights,
glmfamily, Rcontrol, Rlaplace,
betapriorfamily));
prior_m = compute_prior_probs(model,pmodel,p, modelprior);
prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);

logmargy = REAL(getListElement(getListElement(glm_fit, "lpy"),"lpY"))[0];
shrinkage_m = REAL(getListElement(getListElement(glm_fit, "lpy"),
Expand Down Expand Up @@ -140,7 +141,7 @@ SEXP glm_mcmcbas(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
glm_fit = PROTECT(glm_FitModel(X, Y, Rmodel_m, Roffset, Rweights,
glmfamily, Rcontrol, Rlaplace,
betapriorfamily));
prior_m = compute_prior_probs(model,pmodel,p, modelprior);
prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);

logmargy = REAL(getListElement(getListElement(glm_fit, "lpy"),"lpY"))[0];
shrinkage_m = REAL(getListElement(getListElement(glm_fit, "lpy"),
Expand Down Expand Up @@ -225,7 +226,7 @@ SEXP glm_mcmcbas(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
glm_fit = PROTECT(glm_FitModel(X, Y, Rmodel_m, Roffset, Rweights,
glmfamily, Rcontrol, Rlaplace,
betapriorfamily));
prior_m = compute_prior_probs(model,pmodel,p, modelprior);
prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);
logmargy = REAL(getListElement(getListElement(glm_fit, "lpy"),"lpY"))[0];
shrinkage_m = REAL(getListElement(getListElement(glm_fit, "lpy"),
"shrinkage"))[0];
Expand Down
5 changes: 3 additions & 2 deletions src/glm_sampleworep.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ SEXP glm_sampleworep(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
struct Var *vars = (struct Var *) R_alloc(p, sizeof(struct Var)); // Info about the model variables.
probs = REAL(Rprobs);
int n = sortvars(vars, probs, p);
int noInclusionIs1 = no_prior_inclusion_is_1(p, probs);

int *model = ivecalloc(p);
/* fill in the sure things */
Expand Down Expand Up @@ -86,7 +87,7 @@ SEXP glm_sampleworep(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
SEXP glm_fit = PROTECT(glm_FitModel(X, Y, Rmodel_m, Roffset, Rweights,
glmfamily, Rcontrol, Rlaplace,
betapriorfamily));
double prior_m = compute_prior_probs(model,pmodel,p, modelprior);
double prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);
logmargy = REAL(getListElement(getListElement(glm_fit, "lpy"),"lpY"))[0];
shrinkage_m = REAL(getListElement(getListElement(glm_fit, "lpy"),
"shrinkage"))[0];
Expand Down Expand Up @@ -122,7 +123,7 @@ SEXP glm_sampleworep(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights,
glm_fit = PROTECT(glm_FitModel(X, Y, Rmodel_m, Roffset, Rweights,
glmfamily, Rcontrol, Rlaplace,
betapriorfamily));
prior_m = compute_prior_probs(model,pmodel,p, modelprior);
prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);
logmargy = REAL(getListElement(getListElement(glm_fit, "lpy"),"lpY"))[0];
shrinkage_m = REAL(getListElement(getListElement(glm_fit, "lpy"),
"shrinkage"))[0];
Expand Down
3 changes: 2 additions & 1 deletion src/lm_deterministic.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,7 @@ SEXP deterministic(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit,
vars = (struct Var *) R_alloc(p, sizeof(struct Var));
probs = REAL(Rprobs);
n = sortvars(vars, probs, p);
int noInclusionIs1 = no_prior_inclusion_is_1(p, probs);


/* Make space for the models and working variables. */
Expand Down Expand Up @@ -154,7 +155,7 @@ SEXP deterministic(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit,
gexpectations(p, rank_m, nobs, R2_m, alpha, INTEGER(method)[0], RSquareFull,
SSY, &logmarg_m, &shrinkage_m);
REAL(logmarg)[m] = logmarg_m;
REAL(priorprobs)[m] = compute_prior_probs( model, pmodel,p, modelprior);
REAL(priorprobs)[m] = compute_prior_probs( model, pmodel,p, modelprior, noInclusionIs1);
REAL(shrinkage)[m] = shrinkage_m;
UNPROTECT(2);
}
Expand Down
5 changes: 3 additions & 2 deletions src/lm_mcmc.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -317,6 +317,7 @@ extern SEXP mcmc_new(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeld
n = sortvars(vars, probs, p);
for (i =n; i <p; i++) REAL(MCMCprobs)[vars[i].index] = probs[vars[i].index];
for (i =0; i <n; i++) REAL(MCMCprobs)[vars[i].index] = 0.0;
int noInclusionIs1 = no_prior_inclusion_is_1(p, probs);


// fill in the sure things
Expand Down Expand Up @@ -359,7 +360,7 @@ extern SEXP mcmc_new(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeld
&shrinkage_m);


prior_m = compute_prior_probs(model,pmodel,p, modelprior);
prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);
if (prior_m == 0.0) Rprintf("warning initial model has 0 prior probabilty\n");
SetModel2(logmargy, shrinkage_m, prior_m, sampleprobs, logmarg, shrinkage, priorprobs, m);
SetModel(Rcoef_m, Rse_m, Rmodel_m, mse_m, R2_m, beta, se, modelspace, mse, R2, m);
Expand Down Expand Up @@ -405,7 +406,7 @@ extern SEXP mcmc_new(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeld
}

if (newmodel == 1) {
prior_m = compute_prior_probs(model,pmodel,p, modelprior);
prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);
if (prior_m == 0.0) {
MH *= 0.0;
}
Expand Down
7 changes: 4 additions & 3 deletions src/lm_mcmcbas.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,7 @@ SEXP mcmcbas(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeldim, SEXP
for (i =n; i <p; i++) REAL(MCMCprobs)[vars[i].index] = probs[vars[i].index];
for (i =0; i <n; i++) REAL(MCMCprobs)[vars[i].index] = 0.0;
MCMC_probs = REAL(MCMCprobs);
int noInclusionIs1 = no_prior_inclusion_is_1(p, probs);


pigamma = vecalloc(p);
Expand Down Expand Up @@ -229,7 +230,7 @@ SEXP mcmcbas(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeldim, SEXP
REAL(sampleprobs)[m] = 1.0;
REAL(logmarg)[m] = logmargy;
REAL(shrinkage)[m] = shrinkage_m;
prior_m = compute_prior_probs(model,pmodel,p, modelprior);
prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);
REAL(priorprobs)[m] = prior_m;

UNPROTECT(3);
Expand Down Expand Up @@ -313,7 +314,7 @@ SEXP mcmcbas(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeldim, SEXP
memcpy(coefficients, XtYwork, sizeof(double)*pmodel);
cholreg(XtYwork, XtXwork, coefficients, se_m, &mse_m, pmodel, nobs);
if (pmodel > 1) R2_m = 1.0 - (mse_m * (double) ( nobs - pmodel))/SSY;
prior_m = compute_prior_probs(model,pmodel,p, modelprior);
prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);
gexpectations(p, pmodel, nobs, R2_m, alpha, INTEGER(method)[0], RSquareFull, SSY, &logmargy, &shrinkage_m);
postnew = logmargy + log(prior_m);
}
Expand Down Expand Up @@ -515,7 +516,7 @@ SEXP mcmcbas(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeldim, SEXP
gexpectations(p, pmodel, nobs, R2_m, alpha, INTEGER(method)[0], RSquareFull, SSY, &logmargy, &shrinkage_m);
REAL(logmarg)[m] = logmargy;
REAL(shrinkage)[m] = shrinkage_m;
REAL(priorprobs)[m] = compute_prior_probs(model,pmodel,p, modelprior);
REAL(priorprobs)[m] = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);


if (m > 1) {
Expand Down
5 changes: 3 additions & 2 deletions src/lm_sampleworep.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -241,6 +241,7 @@ extern SEXP sampleworep_new(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit,
struct Var *vars = (struct Var *) R_alloc(p, sizeof(struct Var)); // Info about the model variables.
probs = REAL(Rprobs);
int n = sortvars(vars, probs, p);
int noInclusionIs1 = no_prior_inclusion_is_1(p, probs);

SEXP Rse_m = NULL, Rcoef_m = NULL, Rmodel_m = NULL;
RSquareFull = CalculateRSquareFull(XtY, XtX, XtXwork, XtYwork, Rcoef_m, Rse_m, p, nobs, yty, SSY);
Expand Down Expand Up @@ -293,7 +294,7 @@ extern SEXP sampleworep_new(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit,
// gexpectations(p, pmodel, nobs, R2_m, alpha, INTEGER(method)[0], RSquareFull, SSY, &logmargy, &shrinkage_m);

// check should this depend on rank or pmodel?
double prior_m = compute_prior_probs(model,pmodel,p, modelprior);
double prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);



Expand Down Expand Up @@ -349,7 +350,7 @@ extern SEXP sampleworep_new(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit,
// Rprintf("rank %d dim %d\n", rank_m, pmodel);
// gexpectations(p, pmodel, nobs, R2_m, alpha, INTEGER(method)[0], RSquareFull, SSY, &logmargy, &shrinkage_m);

prior_m = compute_prior_probs(model,pmodel,p, modelprior);
prior_m = compute_prior_probs(model,pmodel,p, modelprior, noInclusionIs1);
SetModel2(logmargy, shrinkage_m, prior_m, sampleprobs, logmarg, shrinkage, priorprobs, m);
SetModel_lm(Rcoef_m, Rse_m, Rmodel_m, mse_m, R2_m, beta, se, modelspace, mse, R2,m);
UNPROTECT(3);
Expand Down
17 changes: 16 additions & 1 deletion src/model_probabilities.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,14 +60,29 @@ void compute_margprobs_old(Bit **models, SEXP Rmodelprobs, double *margprobs, in
}
}

double compute_prior_probs(int *model, int modeldim, int p, SEXP modelprior) {
int no_prior_inclusion_is_1(int p, double *probs) {

int noInclusionIs1 = 0;
// loop starts from 1 since the intercept is corrected for in the model prior functions
for (int i = 1; i < p; i++) {
if (probs[i] > (1.0 - DBL_EPSILON)) {
noInclusionIs1++;
}
}
return noInclusionIs1;
}

double compute_prior_probs(int *model, int modeldim, int p, SEXP modelprior, int noInclusionIs1) {
const char *family;
double *hyper_parameters, priorprob = 1.0;


family = CHAR(STRING_ELT(getListElement(modelprior, "family"),0));
hyper_parameters = REAL(getListElement(modelprior,"hyper.parameters"));

// reduce the model space by the number of predictors that are always included
p -= noInclusionIs1;
modeldim -= noInclusionIs1;

if (strcmp(family, "Beta-Binomial") == 0)
priorprob = beta_binomial(modeldim, p, hyper_parameters);
Expand Down
11 changes: 11 additions & 0 deletions tests/testthat/test-model-priors.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,3 +37,14 @@ test_that("Bernoulli hereditary prior", {
modelprior = Bernoulli.heredity(.5, NULL))
)
})

test_that("Always include changes model-prior", {

res <- bas.lm(Y ~ ., data = Hald, modelprior = beta.binomial(1, 1),
include.always = ~ 1 + X1 + X2)
ord <- order(lengths(res$which))
expect_equal(
object = res$priorprobs[ord],
expected = c(0.333333333333333, 0.166666666666667, 0.166666666666667, 0.333333333333333)
)
})