rm(ls = list())
x = 2
rm(ls = list())
rm(list = ls())
library(glmnet)
library(monomvn)
set.seed(2025)
setwd("~/Downloads/diabetes")
# Dados do estudo de diabetes
dados <- read.table('./diabetesdata.txt', header = TRUE)
data_quadmodel <- read.table('./diabetes_data_quadmodel.txt', header = TRUE)
# View(dados)
# View(data_quadmodel)
dim(dados)
dim(data_quadmodel)
names(dados)
###################################
# Organizando os dados para o modelo
n <- nrow(dados)
# Dados do estudo de diabetes
dados <- read.table('./diabetesdata.txt', header = TRUE)
# View(dados)
# View(data_quadmodel)
dim(dados)
dim(data_quadmodel)
names(dados)
###################################
# Organizando os dados para o modelo
n <- nrow(dados)
p <- ncol(dados) - 1
# normalizando os dados
x <- matrix(NA, nrow=n, ncol=p)
for(i in 1:p){
x[,i] <- (dados[,i] - mean(dados[,i]))/(sd(dados[,i])*sqrt(n-1))
}
rm(list = ls())
library(glmnet)
library(monomvn)
set.seed(2025)
# Dados do estudo de diabetes
dados <- read.table('./diabetesdata.txt', header = TRUE)
# View(dados)
# View(data_quadmodel)
dim(dados)
dim(data_quadmodel)
# View(dados)
# View(data_quadmodel)
dim(dados)
names(dados)
###################################
# Organizando os dados para o modelo
n <- nrow(dados)
p <- ncol(dados) - 1
# normalizando os dados
x <- matrix(NA, nrow=n, ncol=p)
for(i in 1:p){
x[,i] <- (dados[,i] - mean(dados[,i]))/(sd(dados[,i])*sqrt(n-1))
}
y <- (dados$Y - mean(dados$Y))/(sd(dados$Y)*sqrt(441))
colnames(x) <- colnames(dados[,1:p])
head(x)
###################################
# ajuste do lasso
lasso_fit <- glmnet::glmnet(x, y, family = 'gaussian', intercept = FALSE,
alpha = 1)
# png(file="diabetes_lasso_path.png",
#      width=600, height=500, res = 100)
plot(lasso_fit, xvar = 'lambda', sign.lambda = 1,
lwd=2, ylab='Coeficientes')
# validação cruzada para escolher lambda
cv_lasso_fit <- glmnet::cv.glmnet(x, y, family = 'gaussian', intercept = FALSE,
alpha = 1, type.measure = 'mse', nfolds = 10)
# png(file="diabetes_lasso_mse_cv.png",
#     width=600, height=500, res = 100)
plot(cv_lasso_fit, sign.lambda = 1)
# coeficientes estimados usando o lambda 'ótimo'
cv_lasso_fit$lambda.1se
id_labda_min <- which(lasso_fit$lambda == cv_lasso_fit$lambda.1se)
beta_hat_lasso <- coef(lasso_fit, s=lasso_fit$lambda[id_labda_min])
beta_hat_lasso
n_mcmc <- 5000 # iteracoes de MCMC
n_thin <- 5 # lag entre os valores mantidos
burnin <- 1000 # numero de valores iniciais descartados
blasso_fit <- blasso(X = x,
y = y,
icept = FALSE,
T = n_mcmc,
thin = n_thin,
verb = 1)
blasso_fit_summary <- summary(blasso_fit)
dim(blasso_fit$beta)
plot(blasso_fit)
beta_mcmc <- blasso_fit$beta[-c(1:burnin),]
lambda2_mcmc <- blasso_fit$lambda2[-c(1:burnin)]
s2_mcmc <- blasso_fit$s2[-c(1:burnin)]
# algumas cadeias de MCMC
dim(beta_mcmc)
# png(file="diabetes_blasso_mcmc.png",
#     width=1000, height=500, res = 100)
par(mfrow=c(1,3), mar = c(4,4,1,1), mgp=c(2.5,1,0))
plot(beta_mcmc[,1], type='l', ylab=expression(beta[1]),
cex.lab=1.5)
plot(beta_mcmc[,2], type='l', ylab=expression(beta[2]),
cex.lab=1.5)
plot(s2_mcmc, type='l', ylab=expression(sigma^2),
cex.lab=1.5)
# histograma da distribuicao a posteriori
# png(file="diabetes_blasso_beta2_hist.png",
#     width=800, height=500, res = 100)
hist(beta_mcmc[,2], prob=TRUE, xlab=expression(beta[2]), ylab='Frequência',
xlim=c(-.3, .3), cex.lab=1.5, main="", col='gray')
dev.off()
# png(file="diabetes_blasso_mcmc.png",
#     width=1000, height=500, res = 100)
# par(mfrow=c(1,3), mar = c(4,4,1,1), mgp=c(2.5,1,0))
plot(beta_mcmc[,1], type='l', ylab=expression(beta[1]),
cex.lab=1.5)
plot(beta_mcmc[,2], type='l', ylab=expression(beta[2]),
cex.lab=1.5)
plot(s2_mcmc, type='l', ylab=expression(sigma^2),
cex.lab=1.5)
# histograma da distribuicao a posteriori
# png(file="diabetes_blasso_beta2_hist.png",
#     width=800, height=500, res = 100)
hist(beta_mcmc[,2], prob=TRUE, xlab=expression(beta[2]), ylab='Frequência',
xlim=c(-.3, .3), cex.lab=1.5, main="", col='gray')
tmp <- seq(-.3, .3, length=1000)
lines(tmp, exp(-abs(tmp)/sd(y))/(2*sd(y)), col='red', lwd=2)
legend("topright", legend=c('posteriori', 'priori'),
fill=c('gray', 'red'), bty='n', cex = 1.5)
# aproximacao da trajetoria do lasso no caso Bayesiano
beta_l1norm = apply(beta_mcmc, 1, function(x) sum(abs(x)))
# png(file="diabetes_blasso_path_mcmc.png",
#     width=800, height=500, res = 100)
par(mar = c(4,4,1,4), mgp=c(2,1,0))
plot(smooth.spline(beta_l1norm, beta_mcmc[,1], spar = 1),
col = 1, lwd = 2, type='l', ylim=c(-.5,.5),
ylab="Coeficientes", xlab=expression("||" * beta * "||"[1]))
for(i in 2:p){
lines(smooth.spline(beta_l1norm, beta_mcmc[,i], spar = 1),
col = i, lwd = 2)
}
legend("right",
inset = c(-0.12, 0),
legend = c(expression(beta[1]),expression(beta[2]),
expression(beta[3]),expression(beta[4]),
expression(beta[5]),expression(beta[6]),
expression(beta[7]),expression(beta[8]),
expression(beta[9]),expression(beta[10])),
fill = c(1:10),
bty='n',
xpd = NA)
# aproximacao da trajetoria do lasso no caso Bayesiano
beta_l1norm = apply(beta_mcmc, 1, function(x) sum(abs(x)))
# png(file="diabetes_blasso_boxplots.png",
#     width=800, height=500, res = 100)
boxplot(beta_mcmc, names=c(expression(beta[1]),expression(beta[2]),
expression(beta[3]),expression(beta[4]),
expression(beta[5]),expression(beta[6]),
expression(beta[7]),expression(beta[8]),
expression(beta[9]),expression(beta[10])),
ylab="Coeficientes")
abline(h=0, lty=2, lwd=2, col='red')
# médias, erro padrão e quantis a posteriori
round(apply(beta_mcmc, 2, mean), 4)
round(apply(beta_mcmc, 2, sd), 4)
round(apply(beta_mcmc, 2, function(x) quantile(x, c(.025, .975))), 4)
## sumário da variância
plot(blasso_fit, burnin=100, which="s2")
blasso_fit_summary$s2
## sumário do parâmetro de penalização
plot(blasso_fit, burnin=200, which="lambda2")
blasso_fit_summary$lambda2
rm(list = ls())
library(glmnet)
set.seed(2025)
# Dados do estudo de diabetes
dados <- read.table('./diabetesdata.txt', header = TRUE)
data_quadmodel <- read.table('./diabetes_data_quadmodel.txt', header = TRUE)
# View(dados)
# View(data_quadmodel)
dim(dados)
dim(data_quadmodel)
names(dados)
rm(list = ls())
library(glmnet)
set.seed(2025)
# Dados do estudo de diabetes
dados <- read.table('./diabetesdata.txt', header = TRUE)
# View(dados)
# View(data_quadmodel)
dim(dados)
names(dados)
###################################
# Organizando os dados para o modelo
n <- nrow(dados)
p <- ncol(dados) - 1
# normalizando os dados
x <- matrix(NA, nrow=n, ncol=p)
for(i in 1:p){
x[,i] <- (dados[,i] - mean(dados[,i]))/(sd(dados[,i])*sqrt(n-1))
}
y <- (dados$Y - mean(dados$Y))/(sd(dados$Y)*sqrt(441))
colnames(x) <- colnames(dados[,1:p])
head(x)
###################################
# ajuste do lasso
lasso_fit <- glmnet::glmnet(x, y, family = 'gaussian', intercept = FALSE,
alpha = 1)
# png(file="diabetes_lasso_path.png",
#      width=600, height=500, res = 100)
plot(lasso_fit, xvar = 'lambda', sign.lambda = 1,
lwd=2, ylab='Coeficientes')
###################################
# ajuste do lasso
lasso_fit <- glmnet::glmnet(x, y, family = 'gaussian', intercept = FALSE,
alpha = 1)
# png(file="diabetes_lasso_path.png",
#      width=600, height=500, res = 100)
plot(lasso_fit, xvar = 'lambda', sign.lambda = 1,
lwd=2, ylab='Coeficientes')
# validação cruzada para escolher lambda
cv_lasso_fit <- glmnet::cv.glmnet(x, y, family = 'gaussian', intercept = FALSE,
alpha = 1, type.measure = 'mse', nfolds = 10)
# png(file="diabetes_lasso_mse_cv.png",
#     width=600, height=500, res = 100)
plot(cv_lasso_fit, sign.lambda = 1)
# coeficientes estimados usando o lambda 'ótimo'
cv_lasso_fit$lambda.1se
beta_hat_lasso <- coef(lasso_fit, s=cv_lasso_fit$lambda.1se)
beta_hat_lasso
B <- 300 # numero de replicas bootstrap
# arrays para salvar as replicas bootstrap
lasso_fit_boot <- list()
lasso_fit_boot$beta_hat_lasso <- matrix(NA, nrow=B, ncol=p)
lasso_fit_boot$lambda1se <- rep(NA, B)
coefpath_boot <- list()
lambda_boot <- list()
for(b in 1:B){
# reamostragem bootstrap
id_boot <- sample(1:n, replace = TRUE)
x_boot <- x[id_boot, ]
y_boot <- y[id_boot]
# modelo estimado bootstrap
cv_lasso_star <- glmnet::cv.glmnet(x_boot, y_boot, family = 'gaussian',
intercept = FALSE, alpha = 1,
type.measure = 'mse', nfolds = 10)
# estimativas boostrap
lambda_star <- cv_lasso_star$lambda.1se
beta_hat_lasso_star <- coef(cv_lasso_star$glmnet.fit, s=lambda_star)
lasso_fit_boot$beta_hat_lasso[b,] <- as.numeric(beta_hat_lasso_star[2:(p+1)])
lasso_fit_boot$lambda1se[b] <- lambda_star
coefpath_boot[[b]] <- cv_lasso_star$glmnet.fit$beta
lambda_boot[[b]] <- cv_lasso_star$glmnet.fit$lambda
# progresso
cat(paste(floor(100*b/B),'%\n'))
}
# png(file="diabetes_boot_beta1_beta2_lambda.png",
#     width=1000, height=500, res = 100)
# par(mfrow=c(1,3), mar = c(4,4,1,1), mgp=c(2.5,1,0))
# beta 1
lmin <- min(lasso_fit_boot$beta_hat_lasso[,1])
lmax <- max(lasso_fit_boot$beta_hat_lasso[,1])
hist(lasso_fit_boot$beta_hat_lasso[,1], prob=TRUE,
xlab=expression(beta[1]), ylab='Frequência',
xlim=c(lmin-.01, lmax+.01), cex.lab=1.5, main="", col='gray')
# beta 2
lmin <- min(lasso_fit_boot$beta_hat_lasso[,2])
lmax <- max(lasso_fit_boot$beta_hat_lasso[,2])
hist(lasso_fit_boot$beta_hat_lasso[,2], prob=TRUE,
xlab=expression(beta[2]), ylab='Frequência',
xlim=c(lmin-.01, lmax+.01), cex.lab=1.5, main="", col='gray')
# lambda
lmin <- min(lasso_fit_boot$lambda1se)
lmax <- max(lasso_fit_boot$lambda1se)
hist(lasso_fit_boot$lambda1se, prob=TRUE,
xlab=expression(lambda), ylab='Frequência',
xlim=c(lmin-.001, lmax+.001), cex.lab=1.5, main="", col='gray')
# png(file="diabetes_boot_lassopath.png",
#     width=1000, height=500, res = 100)
# par(mfrow=c(1,2), mar = c(4,4,1,1), mgp=c(2.5,1,0))
# beta 1
plot(-log(lasso_fit$lambda), lasso_fit$beta[1,], type='n',
xlab=expression(-log(lambda)), ylab=expression(beta[1]), cex.lab=1.5,
ylim=c(-.1, .1))
for(b in 1:B){
lines(-log(lambda_boot[[b]]), coefpath_boot[[b]][1,], col='LightGray')
}
lines(-log(lasso_fit$lambda), lasso_fit$beta[1,], col='Red', lwd=2)
abline(v=-log(cv_lasso_fit$lambda.1se), lty=2, lwd=2, col=1)
# beta 2
plot(-log(lasso_fit$lambda), lasso_fit$beta[2,], type='n',
xlab=expression(-log(lambda)), ylab=expression(beta[2]), cex.lab=1.5,
ylim=c(-.25, .1))
for(b in 1:B){
lines(-log(lambda_boot[[b]]), coefpath_boot[[b]][2,], col='LightGray')
}
lines(-log(lasso_fit$lambda), lasso_fit$beta[2,], col='Red', lwd=2)
abline(v=-log(cv_lasso_fit$lambda.1se), lty=2, lwd=2, col=1)
# png(file="diabetes_boot_boxplots.png",
#     width=800, height=500, res = 100)
boxplot(lasso_fit_boot$beta_hat_lasso,
names=c(expression(beta[1]),expression(beta[2]),
expression(beta[3]),expression(beta[4]),
expression(beta[5]),expression(beta[6]),
expression(beta[7]),expression(beta[8]),
expression(beta[9]),expression(beta[10])),
ylab="Coeficientes", cex.lab=1.5)
abline(h=0, lty=2, lwd=2, col='red')
# png(file="diabetes_boot_probzero.png",
#     width=800, height=500, res = 100)
barplot(colMeans(lasso_fit_boot$beta_hat_lasso==0),
ylab = 'Probabilidade de ser zero',
names=c(expression(beta[1]),expression(beta[2]),
expression(beta[3]),expression(beta[4]),
expression(beta[5]),expression(beta[6]),
expression(beta[7]),expression(beta[8]),
expression(beta[9]),expression(beta[10])),
cex.lab=1.5)
# erros padrão bootstrap
round(apply(lasso_fit_boot$beta_hat_lasso, 2, sd), 4)
# intervalos de confiança bootstrap percentil
round(apply(lasso_fit_boot$beta_hat_lasso, 2,
function(x) quantile(x, c(.025, .975))), 4)
rm(list = ls())
library(glmnet)
source("lasso_inference.R")
set.seed(2025)
# Dados do estudo de diabetes
dados <- read.table('./diabetesdata.txt', header = TRUE)
data_quadmodel <- read.table('./diabetes_data_quadmodel.txt', header = TRUE)
# View(dados)
# View(data_quadmodel)
dim(dados)
dim(data_quadmodel)
names(dados)
# dados com variaveis do modelo quadratico normalizadas
x <- as.matrix(data_quadmodel)
colnames(x) <- colnames(data_quadmodel)
colnames(x)[1:15]
# normalizando os dados da variavel resposta
y <- (dados$Y - mean(dados$Y))/(sd(dados$Y)*sqrt(441))
# dimensao dos dados
n <- nrow(x)
p <- ncol(x)
###################################
# ajuste do lasso
lasso_fit <- glmnet::glmnet(x, y, family = 'gaussian', intercept = FALSE,
alpha = 1)
# png(file="diabetes_quad_lasso_path.png",
#      width=600, height=500, res = 100)
plot(lasso_fit, xvar = 'lambda', sign.lambda = 1,
lwd=2, ylab='Coeficientes')
# validação cruzada para escolher lambda
cv_lasso_fit <- glmnet::cv.glmnet(x, y, family = 'gaussian', intercept = FALSE,
alpha = 1, type.measure = 'mse', nfolds = 10)
# png(file="diabetes_quad_lasso_mse_cv.png",
#     width=600, height=500, res = 100)
plot(cv_lasso_fit, sign.lambda = 1)
# coeficientes estimados usando o lambda 'ótimo'
cv_lasso_fit$lambda.min
beta_hat_lasso <- coef(lasso_fit, s=cv_lasso_fit$lambda.min)
sum(beta_hat_lasso!=0)
# png(file="diabetes_quad_lasso_coefs.png",
#     width=600, height=500, res = 100)
plot(beta_hat_lasso[2:(p+1)], pch=' ', cex.lab=1.4,
cex.axis=1.5, cex=1, lwd = 1.5, ylab='Estimativa', xlab='índice')
grid()
abline(h=0, lwd=0.5, lty=2, col='red')
points(beta_hat_lasso[2:(p+1)], col="black", pch=20)
# dez maiores estimativas pelo lasso
id_top10 <- order(abs(beta_hat_lasso), decreasing = TRUE)[1:10]
round(beta_hat_lasso[id_top10], 4)
colnames(x)[id_top10]
devlasso_fit <- SSLasso(x,y,verbose = TRUE, alpha = .05, intercept = FALSE,
lambda = cv_lasso_fit$lambda.1se)
lmin <- max(devlasso_fit$up.lim)
lmax <- min(devlasso_fit$low.lim)
# png("diabetes_quad_devlasso_coefs.png",
#     width=600, height=500, res = 100)
plot(devlasso_fit$coef, ylim=c(lmin, lmax), pch=20, cex.lab=1.4,
cex.axis=1.5, cex=1, lwd = 1.5, ylab='Estimativa', xlab='índice')
grid()
for(j in 1:p){
segments(j,devlasso_fit$low.lim[j],j,devlasso_fit$up.lim[j],
col="orange",lwd=2)
}
abline(h=0, lwd=0.5, lty=2, col='red')
points(devlasso_fit$unb.coef,col="blue", pch=20)
points(devlasso_fit$coef, col="black", pch=20)
legend("bottomright",c("Lasso","Lasso deviesado","IC de 95%"),
col=c("black","blue","orange"),pch=c(20,20,NA_integer_),
lty = c(0,0,1), bty='n', cex=1.1,lwd=2)
# indices dos ICs que contem zero mas cuja estimativa lasso é diferente de zero
lasso_vs_devlasso <- (devlasso_fit$coef!=0)*(devlasso_fit$up.lim>0)*
(devlasso_fit$low.lim<0)
id_dif <- which(lasso_vs_devlasso == 1)
id_dif
# png("diabetes_quad_lassodif0_coefs_IC_comzero.png",
#     width=600, height=500, res = 100)
plot(devlasso_fit$coef, ylim=c(lmin, lmax), pch=' ', cex.lab=1.4,
cex.axis=1.5, cex=1, lwd = 1.5, ylab='Estimativa', xlab='índice')
grid()
for(j in id_dif){
segments(j,devlasso_fit$low.lim[j],j,devlasso_fit$up.lim[j],
col="orange",lwd=2)
}
abline(h=0, lwd=0.5, lty=2, col='red')
points(id_dif, devlasso_fit$unb.coef[id_dif],col="blue", pch=20)
points(id_dif, devlasso_fit$coef[id_dif], col="black", pch=20)
legend("bottomright",c("Lasso","Lasso deviesado","IC de 95%"),
col=c("black","blue","orange"),pch=c(20,20,NA_integer_),
lty = c(0,0,1), bty='n', cex=1.1,lwd=2)
# png("diabetes_quad_devlasso_pvalor.png",
#     width=600, height=500, res = 100)
plot(devlasso_fit$pvals, type='h', cex.lab=1.4,
cex.axis=1.5, cex=1, lwd = 1.5, ylab='Valor-p', xlab='índice')
abline(h=0.05, lty=2, col='red', lwd=2)
# dev.off()
# quais covariaveis ficaram com valor-p abaixo de 5%
which(devlasso_fit$pvals < 0.05)
# comparação do top-10 do lasso com valores pelo lasso deviesado
id_top10 <- order(abs(beta_hat_lasso), decreasing = TRUE)[1:10]
round(beta_hat_lasso[id_top10], 4)
round(devlasso_fit$unb.coef[id_top10], 4)
round(devlasso_fit$low.lim[id_top10], 4)
round(devlasso_fit$up.lim[id_top10], 4)
colnames(x)[id_top10]
