CÓDIGO CIGARRILLOS

x <- cigarrillos$age
y <- cigarrillos$income
df2 = data.frame(x, y)

vars <- length(cigarrillos)-1
vars

n <- length(cigarrillos$age) / vars
n

dfree <- n - vars - 1
dfree

Sxy <- sum(x*y) - n * mean(x) * mean(y) 
Sxx <- sum(x^2) - n* (mean(x))^2
Syy <- sum(y^2) - n* (mean(y))^2

cat( " Sxy = ", Sxy,"\n Sxx = ", Sxx, "\n Syy =",Syy)


r <- Sxy / sqrt(Sxx * Syy)
cat( "Coeficiente r = ", r)


t.value <- (r * sqrt(n-2)) / sqrt(1-r^2)

Valor-p para r 

p.value <- pt(q = t.value, df = n-2, lower.tail = FALSE)

cat("Estadístico T=", round(t.value,3), 
    "\nValor P = ", round(p.value, 5))



Regresión lineal 

Parámetros de regresión

b1 <- Sxy / Sxx
b0 <- mean(y) - b1 * mean(x)

cat(" b0 = ", b0, "\n b1 = ", b1)

ANOVA Y SUMA DE CUADRADOS

SCM <- b1*Sxy
SCE <- Syy - b1*Sxy


CMM <- SCM / vars
CME <- SCE / (n-2)


fstat <- CMM / CME
fpvalue <- pf(fstat , 1, n-2, lower.tail = FALSE)



cat(
  "Fuente de |", "\t Grados de ", "\t Sumas de", "\t Cuadrados",
  "\nVariación |", "\t Libertad  ", "\t Cuadrados", "\t Medios",
  
  "\nModelo    | \t", vars, "\t \t" , SCM,"\t " , CMM,
  "\nResidual  | \t",dfree, "\t \t" , SCE,"\t " ,CME,
  "\nTotal     | \t", (vars + dfree), "\t \t", Syy)

cat(" F = ", fstat, ", Valor P = ", fpvalue)

b1.err <- sqrt(CME / Sxx)
b0.err <- sqrt(CME * sum(x^2) / (n*Sxx))

b0.t <- (b0 - 0) / b0.err
b1.t <- (b1 - 0) / b1.err

b0.p <- 2 * pt(b0.t, df = n-2, lower.tail = FALSE)
b1.p <- 2 * pt(b1.t, df = n-2, lower.tail = FALSE)


cat(" T_b0= ",b0.t, "Valor P = ", b0.p,
    "\n T_b1= ", b1.t, "Valor P = ", b1.p)



r2 <- (Syy- SCE)/Syy



Análisis de residuales

fitted <- x * b1 + b0
resd <- y - fitted
min.res <- round(min(resd), 3)
max.res <- round(max(resd), 3)
q1.q3 <- quantile(resd, probs = c(.25, .75))
med <- round(median(resd), 3)



coeffs <- data.frame(
  cbind(c(b0, b1), c(b0.err, b1.err), c(b0.t, b1.t),c(b0.p, b1.p))
)
colnames(coeffs) <- c('Estimate', 'Std. Error', 't value', 'Pr(>|t|)')

rownames(coeffs) <- c('Intercept', 'x1')

rsquare <- paste('Multiple R-squared: ', round(r2, 4))
f.stat <- paste(
  'Estadistico F: ', round(fstat, 2), ' con 1 y ', n - 2,
  ' G.L., Valor P: ', format(fpvalue, digits = 3, scientific = TRUE)
)
residual <- data.frame(
  cbind(min.res, round(q1.q3[1], 3), med, round(q1.q3[2], 3), max.res)
)
colnames(residual) <- c('Min', 'Q1', 'Median', 'Q3', 'Max')

resdi <- paste('Error residual estándar: ', round(CME, 2), 
               ' con ', n-2, ' grados de libertad')

### Lista con los resúmenes
regresion <- list('Residuales'=residual, 
                  'Coeficientes'=coeffs, 
                  'ErrorResidual'= resdi, 
                  'Determinacion'=rsquare, 
                  'EstadisticoF'=f.stat)

regresion
regresion$Residuales




Predecir valores y bandas de confianza y predicción

predicciones_conf <- predict(modelo, interval = "confidence", level = 0.95)
predicciones_pred <- predict(modelo, interval = "prediction", level = 0.95)

Agregar las predicciones al dataframe

df$prediccion <- predicciones_conf[, 1]
df$lim_inf_conf <- predicciones_conf[, 2]
df$lim_sup_conf <- predicciones_conf[, 3]
df$lim_inf_pred <- predicciones_pred[, 2]
df$lim_sup_pred <- predicciones_pred[, 3]

# Crear el gráfico de dispersión y la línea de regresión con bandas de confianza y predicción
grafico <- ggplot(df, aes(x = x, y = y)) +
  geom_vline(xintercept = 0, linetype = "solid", size = 1.1, color = "grey70") +  # línea del eje x
  geom_hline(yintercept = 0, linetype = "solid", size = 1.1, color = "grey70") +  # línea del eje y
  geom_point(color = "blue") +  # puntos de dispersión
  geom_ribbon(aes(ymin = lim_inf_conf, ymax = lim_sup_conf), fill = "steelblue", alpha = 0.4) +  # banda de confianza
  geom_ribbon(aes(ymin = lim_inf_pred, ymax = lim_sup_pred), fill = "lightgreen", alpha = 0.2) +  # banda de predicción
  geom_line(aes(y = prediccion), color = "green4",  size=1.1, linetype = "solid", alpha = 0.7) +  # línea de regresión
  labs(title = "Regresión lineal con bandas de confianza y de predicciones",
       subtitle = "Velocidad vs Distancia de frenado",# con Bandas de Confianza ",
       x = "Variable Independiente",
       y = "Variable Dependiente")+
  xlab("Speed") +  # etiqueta eje x
  ylab("Distance") +  # etiqueta eje y
  theme_bw()
print(grafico) # Mostrar el gráfico
predicciones_conf
predicciones_pred
head(df)