Post date: Oct 15, 2019 5:14:43 PM
I tried Bayesian ridge and lasso regression as alternatives to the model selection/averaging to model expected fitness. These are two additional methods of penalized/shrinkage inducing regression that can help stabilize things (avoid over fitting) when the number of parameters is somewhat large relative to the number of observations (but not necessarily when p >> n). Ridge regression tends to shrink all coefficients in a more uniform way, wheras lasso shrinks the weakest effects most strongly. This is a nice post that I based my code and approach on: http://haines-lab.com/post/on-the-equivalency-between-the-lasso-ridge-regression-and-specific-bayesian-priors/.
My analysis used rstan (Version 2.19.2, GitRev: 2e1f913d3ca3) and R 3.5.2 (on my computer).
This uses a subset of code from models.R that has been preserved in /uufs/chpc.utah.edu/common/home/u6000989/projects/timema_genarch/color_exp_2019_gbs/Gemma_mel/ridgeLassoTchum/ and is pasted below.
load("models.rdat")
## bayesian ridge
fit_AC<-read.table("exfit_AC.txt",header=TRUE)
fit_MM<-read.table("exfit_MM.txt",header=TRUE)
library(rstan)
modelStringRidge<-"data{
int N_train; // # training observations
int N_pred; // # predictor variables
vector[N_train] y_train; // training outcomes
matrix[N_train, N_pred] X_train; // training data
}
parameters{
real<lower=0> sigma; // error SD
real<lower=0> sigma_B; // hierarchical SD across betas
vector[N_pred] beta; // regression beta weights
}
model{
// group-level (hierarchical) SD across betas
sigma_B ~ cauchy(0, 1);
// model error SD
sigma ~ normal(0, 1);
// beta prior (provides ridge regularization)
beta ~ normal(0, sigma_B);
// model likelihood
y_train ~ normal(X_train*beta, sigma);
}
"
stanDsoRidge<-stan_model(model_code=modelStringRidge)
modelStringLasso<-"data{
int N_train; // # training observations
int N_pred; // # predictor variables
vector[N_train] y_train; // training outcomes
matrix[N_train, N_pred] X_train; // training data
}
parameters{
real<lower=0> sigma; // error SD
real<lower=0> sigma_B; // hierarchical SD across betas
vector[N_pred] beta; // regression beta weights
}
model{
// group-level (hierarchical) SD across betas
sigma_B ~ cauchy(0, 1);
// model error SD
sigma ~ normal(0, 1);
// beta prior (provides lasso regularization)
beta ~ double_exponential(0, sigma_B);
// model likelihood
y_train ~ normal(X_train*beta, sigma);
}
"
stanDsoLasso<-stan_model(model_code=modelStringLasso)
V<-31 ## variables
form<-p_AC[,2] ~ .^5
mod<-model.matrix(form,dat=as.data.frame(c_AC_pip[,-6]))
mod<-cbind(mod[,-1])
X<-data.frame(as.numeric(fit_AC$blk==1),as.numeric(fit_AC$blk==2),as.numeric(fit_AC$blk==3),
mod)
nas<-which(is.na(fit_AC$exfit)==TRUE)
kp<-1:219
kp<-kp[-nas]
sl_AC<-list(N_train=216,N_pred=V+3,
y_train=fit_AC$exfit[kp], X_train=X[kp,])
fit_ridge_AC <- sampling(stanDsoRidge, sl_AC, iter = 2000,
warmup = 500, chains = 3, cores = 3)
fit_lasso_AC <- sampling(stanDsoLasso, sl_AC, iter = 2000,
warmup = 500, chains = 3, cores = 3)
est_ridge_AC<-rstan::extract(fit_ridge_AC)
est_lasso_AC<-rstan::extract(fit_lasso_AC)
beta_lasso_AC<-apply(est_lasso_AC[[3]],2,mean)[-c(1:3)]
beta_ridge_AC<-apply(est_ridge_AC[[3]],2,mean)[-c(1:3)]
V<-31 ## variables
form<-p_MM[,2] ~ .^5
mod<-model.matrix(form,dat=as.data.frame(c_MM_pip[,-6]))
mod<-cbind(mod[,-1])
X<-data.frame(as.numeric(fit_MM$blk==1),as.numeric(fit_MM$blk==2),as.numeric(fit_MM$blk==3),
mod)
nas<-which(is.na(fit_MM$exfit)==TRUE)
kp<-1:216
kp<-kp[-nas]
sl_MM<-list(N_train=213,N_pred=V+3,
y_train=fit_MM$exfit[kp], X_train=X[kp,])
fit_ridge_MM <- sampling(stanDsoRidge, sl_MM, iter = 2000,
warmup = 500, chains = 3, cores = 3)
fit_lasso_MM <- sampling(stanDsoLasso, sl_MM, iter = 2000,
warmup = 500, chains = 3, cores = 3)
est_ridge_MM<-rstan::extract(fit_ridge_MM)
est_lasso_MM<-rstan::extract(fit_lasso_MM)
beta_lasso_MM<-apply(est_lasso_MM[[3]],2,mean)[-c(1:3)]
beta_ridge_MM<-apply(est_ridge_MM[[3]],2,mean)[-c(1:3)]
cs<-rep(c("cornsilk3","coral3","gold1","dodgerblue4","black"),c(5,10,10,5,1))
pdf("ColorExFitLassoEffs5snpAllEpi.pdf",width=7,height=9)
par(mfrow=c(2,1))
par(mar=c(5.75,4.5,3,0.5))
barplot(beta_lasso_AC,ylim=c(-.2,.3),names.arg=paIds,col=cs,xlab="Coefficient",ylab="Effect",cex.lab=1.4,las=2,cex.names=.7)
title(main="(A) Expected fitness on A/C",cex.main=1.5)
barplot(beta_lasso_MM,ylim=c(-.2,.3),names.arg=paIds,col=cs,xlab="Coefficient",ylab="Effect",cex.lab=1.4,las=2,cex.names=.7)
title(main="(B) Expected fitness on MM",cex.main=1.5)
dev.off()
pdf("ColorExFitRidgeEffs5snpAllEpi.pdf",width=7,height=9)
par(mfrow=c(2,1))
par(mar=c(5.75,4.5,3,0.5))
barplot(beta_ridge_AC,ylim=c(-.2,.3),names.arg=paIds,col=cs,xlab="Coefficient",ylab="Effect",cex.lab=1.4,las=2,cex.names=.7)
title(main="(A) Expected fitness on A/C",cex.main=1.5)
barplot(beta_ridge_MM,ylim=c(-.2,.3),names.arg=paIds,col=cs,xlab="Coefficient",ylab="Effect",cex.lab=1.4,las=2,cex.names=.7)
title(main="(B) Expected fitness on MM",cex.main=1.5)
dev.off()