angle=90, code=3, length=0.05, col='blue', lwd=1.5)
lambda_cv <- vlambdas[which.min(eqm_vlambdas)]
abline(v=lambda_cv, lty=2, lwd = 3, col='green')
# escolhendo a penalização para o lasso via validação cruzada
cv_lasso_fit <- cv.glmnet(Xpad, ypad, type.measure = 'mse', nfolds = K)
cv_lasso_fit$lambda.min
lambda_cv
png(file="eqm_cv_lambdas.png",
width=600, height=500, res = 100)
plot(vlambdas, eqm_vlambdas, ylim = c(min(ic_low), max(ic_upp)),
ylab='EQM', xlab='lambdas')
arrows(vlambdas,ic_low, vlambdas, ic_upp,
angle=90, code=3, length=0.05, col='blue', lwd=1.5)
dev.off()
# lambda que minimiza o EQM de previsao
lambda_cv <- vlambdas[which.min(eqm_vlambdas)]
lambda_cv
png(file="eqm_cv_lambdas.png",
width=600, height=500, res = 100)
plot(vlambdas, eqm_vlambdas, ylim = c(min(ic_low), max(ic_upp)),
ylab='EQM', xlab='lambdas')
arrows(vlambdas,ic_low, vlambdas, ic_upp,
angle=90, code=3, length=0.05, col='blue', lwd=1.5)
abline(v=lambda_cv, lty=2, lwd = 3, col='green')
dev.off()
max(t(Xpad)%*%ypad)/n
limiar_suave = function(x, lambda){
sinal_x <- ifelse(x>0, 1, -1)
return( sinal_x*max(0, abs(x) - lambda) )
}
x <- seq(-3, 3, length=1000)
limiar_suave = function(x, lambda){
sinal_x <- 1*(x>0)
return( sinal_x*max(0, abs(x) - lambda) )
}
x <- seq(-3, 3, length=10)
limiar_suave(x, 2)
limiar_suave = function(x, lambda){
sinal_x <- 2*(x>0) - 1
return( sinal_x*max(0, abs(x) - lambda) )
}
x <- seq(-3, 3, length=10)
limiar_suave(x, 2)
x
lambda = 2
sinal_x*max(0, abs(x) - lambda)
sinal_x <- 2*(x>0) - 1
sinal_x
abs(x)
max(0, abs(x) - lambda)
max(rep(0, length(x)), abs(x) - lambda)
sapply(x, function(y) max(y, abs(y)-lambda))
sapply(x, function(y) max(0, abs(y)-lambda))
limiar_suave = function(x, lambda){
sinal_x <- 2*(x>0) - 1
return( sinal_x*sapply(x, function(y) max(0, abs(y)-lambda)) )
}
lambda = 2
x <- seq(-3, 3, length=10)
limiar_suave(x, 2)
lambda = 1
x <- seq(-3, 3, length=10)
limiar_suave(x, 1)
plot(x, limiar_suave(x, 1), lty=2)
x <- seq(-3, 3, length=10000)
plot(x, limiar_suave(x, 1), lty=2)
x <- seq(-3, 3, length=1000)
plot(x, limiar_suave(x, 1), lty=2)
plot(x, limiar_suave(x, 1), type='l', lty=2)
abline(a=0)
abline(b=0)
abline(a=0,b=0)
abline(a=0,b=1)
plot(x, limiar_suave(x, 1), type='l', lty=2)
abline(a=0,b=1)
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2, col='blue')
abline(a=0,b=1, lwd=2)
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2, col='blue',
ylab = expression(cal(S)_lambda))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2, col='blue',
ylab = expression(cal(S)[lambda]))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2, col='blue',
ylab = expression(S[lambda](x)))
abline(a=0,b=1, lwd=2)
arrows(1.5,0.5, 1.5, 1.5,
angle=90, code=3, length=0.05, col='blue', lwd=1, lty=3)
text(1.7, .7, expression(lambda), col='blue')
text(1.7, 1.7, expression(lambda), col='blue')
text(1.7, 1.5, expression(lambda), col='blue')
text(1.7, 1.3, expression(lambda), col='blue')
text(1.7, 1.2, expression(lambda), col='blue')
text(1.7, 1.2, expression(lambda), col='blue', cex=2)
text(1.7, 1.2, expression(lambda), col='blue', cex=1.5)
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2, col='blue',
ylab = expression(S[lambda](x)))
abline(a=0,b=1, lwd=2)
arrows(1.5,0.5, 1.5, 1.5,
angle=90, code=3, length=0.05, col='blue', lwd=1, lty=3)
text(1.7, 1.2, expression(lambda), col='blue', cex=1.5)
abline(v=1, lwd=1, lty=3)
abline(v=-11, lwd=1, lty=3)
abline(v=-11, lwd=1, lty=3)
abline(v=-1, lwd=1, lty=3)
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)))
abline(a=0,b=1, lwd=2, col='blue')
arrows(1.5,0.5, 1.5, 1.5,
angle=90, code=3, length=0.05, lwd=1, lty=3)
text(1.7, 1.2, expression(lambda), cex=1.5)
abline(v=1, lwd=1, lty=3)
abline(v=-1, lwd=1, lty=3)
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)), cex.lab=1.5, cex.axis=1.5)
par(mgp = c(4, 1, 0), mar=c(5,5.5,2,0))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)), cex.lab=1.5, cex.axis=1.5)
par(mgp = c(4, 1, 1), mar=c(5,5.5,2,0))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)), cex.lab=1.5, cex.axis=1.5)
par(mar=c(5,5.5,2,0))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)), cex.lab=1.5, cex.axis=1.5)
par(mgp = c(4, 0, 1), mar=c(5,5.5,2,0))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)), cex.lab=1.5, cex.axis=1.5)
par(mgp = c(4, 1, 0), mar=c(5,5.5,2,0))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)), cex.lab=1.5, cex.axis=1.5)
par(mgp = c(4, 1, 0), mar=c(5,5.5,2,1))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)), cex.lab=1.5, cex.axis=1.5)
par(mgp = c(4, 1, 0), mar=c(5,5.5,1,1))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)), cex.lab=1.5, cex.axis=1.5)
par(mgp = c(2, 1, 0), mar=c(5,5.5,1,1))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)), cex.lab=1.5, cex.axis=1.5)
par(mgp = c(3, 1, 0), mar=c(5,5.5,1,1))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)), cex.lab=1.5, cex.axis=1.5)
abline(a=0,b=1, lwd=2, col='blue')
arrows(1.5,0.5, 1.5, 1.5,
angle=90, code=3, length=0.05, lwd=1, lty=3)
text(1.7, 1.2, expression(lambda), cex=1.5)
abline(v=1, lwd=1, lty=3)
abline(v=-1, lwd=1, lty=3)
png(file="limiar_suave.png",
width=600, height=500, res = 100)
par(mgp = c(3, 1, 0), mar=c(5,5.5,1,1))
plot(x, limiar_suave(x, 1), type='l', lty=2, lwd=2,
ylab = expression(S[lambda](x)), cex.lab=1.5, cex.axis=1.5)
abline(a=0,b=1, lwd=2, col='blue')
arrows(1.5,0.5, 1.5, 1.5,
angle=90, code=3, length=0.05, lwd=1, lty=3)
text(1.7, 1.2, expression(lambda), cex=1.5)
abline(v=1, lwd=1, lty=3)
abline(v=-1, lwd=1, lty=3)
dev.off()
library(glmnet)
library(plotmo)
dados <- read.table('./crime.txt')
head(dados)
dim(dados)
# X1 = total overall reported crime rate per 1 million residents
# X2 = reported violent crime rate per 100,000 residents
# X3 = annual police funding in $/resident
# X4 = % of people 25 years+ with 4 yrs. of high school
# X5 = % of 16 to 19 year-olds not in highschool and not highschool graduates.
# X6 = % of 18 to 24 year-olds in college
# X7 = % of people 25 years+ with at least 4 years of college
colnames(dados) <- c('taxa_crime', 'taxa_crime_violento', 'policia',
'em', 'fora_em', 'facul',
'facul4')
summary(dados)
# variáveis
y <- dados$taxa_crime
X <- dados[,c(3,4,5,6,7)]
n <- dim(X)[1]
p <- dim(X)[2]
nrow(X)
ncol(X)
p <- ncol(X)
# variáveis padronizadas
ypad <- y - mean(y)
Xpad <- matrix(0, nrow=n, ncol=p)
for (j in 1:p) {
Xpad[,j] <- (X[,j] - mean(X[,j]))/sd((X[,j] - mean(X[,j])))
}
colnames(Xpad) <- colnames(X)
# escolhendo a penalização para o lasso via validação cruzada
cv_lasso_fit <- cv.glmnet(Xpad, ypad, type.measure = 'mse', nfolds = 10)
cv_lasso_fit$lambda.min
# png(file="crime_mse_cv_lambdas.png",
#    width=600, height=500, res = 100)
plot(cv_lasso_fit, sign.lambda = 1)
png(file="crime_mse_cv_lambdas.png",
width=600, height=500, res = 100)
plot(cv_lasso_fit, sign.lambda = 1)
dev.off()
# Ajuste do lasso
lasso_fit <- glmnet(Xpad, ypad, family = 'gaussian', alpha = 1,
intercept = FALSE)
coef(lasso_fit, s = cv_lasso_fit$lambda.min)
# png(file="trajetoria_lasso_estim.png",
#     width=600, height=600, res = 100)
plotmo::plot_glmnet(lasso_fit, ylab='coeficientes', lwd=2,
main='lasso', cex.main=2)
# calculando erro padrão para as estimativas do lasso usando bootstrap
B <- 500
boot_coef <- matrix(NA, nrow = B, ncol = p)
for (b in 1:B) {
idx <- sample(1:n, replace = TRUE)
fit_b <- glmnet(Xpad[idx, ], ypad[idx], family = 'gaussian', alpha = 1,
intercept = FALSE)
# [-1] remove o intercepto
boot_coef[b, ] <- as.vector(coef(fit_b, s = cv_lasso_fit$lambda.min)[-1])
}
# erros padrão bootstrap
boot_ep <- apply(boot_coef, 2, sd)
boot_ep
# Ajuste de uma regressão ridge
ridge_fit <- glmnet(Xpad, ypad, family = 'gaussian', alpha = 0,
intercept = FALSE)
# png(file="trajetoria_ridge_estim.png",
#     width=600, height=600, res = 100)
plotmo::plot_glmnet(ridge_fit, ylab='coeficientes', lwd=2,
main='Ridge', cex.main=2)
# Ajuste de mínimos quadrados ordinários
mqo_fit <- lm(ypad ~ -1 + Xpad)
sumario_mqo_fit <- summary(mqo_fit)
# Ajuste de mínimos quadrados ordinários usando somente as
# variáveis selecionadas com o lasso
mqo_fit <- lm(ypad ~ -1 + Xpad)
var_selec_lasso <- which(coef(lasso_fit, s = cv_lasso_fit$lambda.min)[-1] != 0)
mqo_fit_relax_lasso <- lm(ypad ~ -1 + Xpad[,var_selec_lasso])
summary(mqo_fit_relax_lasso)
sumario_mqo_fit_relax_lasso <- summary(mqo_fit_relax_lasso)
round(
cbind(sumario_mqo_fit$coefficients[,1:2],
coef(lasso_fit, s = cv_lasso_fit$lambda.min)[-1], boot_ep)
,2)
round(
sumario_mqo_fit_relax_lasso$coefficients[,1:2]
,2)
boot_ep
round(
cbind(sumario_mqo_fit$coefficients[,1:2],
coef(lasso_fit, s = cv_lasso_fit$lambda.min)[-1], boot_ep)
,2)
round(
sumario_mqo_fit_relax_lasso$coefficients[,1:2]
,2)
round(
cbind(sumario_mqo_fit$coefficients[,1:2],
coef(lasso_fit, s = cv_lasso_fit$lambda.min)[-1], boot_ep)
,2)
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)
setwd("C:/Users/rodneyfonseca/Dropbox/Conferences/2026/EBEST/minicurso/diabetes")
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)
###################################
# 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
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)
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]
# 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]
# 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]
