PI2
Sara Gallego Villada y Nicolás Alberto Moreno Reyes
2024-02-23
IMPORTAR DATOS
WORKSPACE
load("C:/Users/sara-/OneDrive - Universidad EAFIT/2024-1/PI2/datos_actualziados.RData")VARIABLES IMPORTANTES
SE FILTRAN PARA DEJAR LAS COMUNAS CON LAS QUE VAMOS A TRABAJAR
letal<-dataframe_prioridad_sin_nulos
comunas_reales_2<- c('Vitacura', 'Las Condes', 'Providencia', 'Lo Barnechea', 'La Reina',
'Nunoa', 'San Miguel', 'Macul', 'Quilicura', 'Santiago', 'Huechuraba',
'Maipu', 'La Florida', 'Penalolen', 'Puente Alto', 'Cerrillos',
'Pudahuel', 'La Cisterna', 'Renca', 'Pedro Aguirre Cerda',
'Estacion Central', 'La Granja', 'Quinta Normal', 'Independencia',
'Recoleta', 'El Bosque', 'San Bernardo', 'San Joaquin', 'Lo Prado',
'Conchali', 'San Ramon', 'Cerro Navia', 'Lo Espejo', 'La Pintana')
# Filtrar el DataFrame para mantener solo las filas con valores permitidos en la primera columna
datos_prioridad_comunas_reales <- subset(dataframe_prioridad_comunas, dataframe_prioridad_comunas$Comuna %in% comunas_reales)
# Filtrar el DataFrame para solo dejar las comunas con valores en hacinamiento
con_comunas <- subset(regiones_y_comunas, regiones_y_comunas$comuna_nom %in% comunas_reales)
hacinamieto_comuna <- con_comunas$hacinamiento
#Regresión de datos funcionales-<Parte matemática
comunas_con_hacinamiento <- con_comunas$comuna_nom
#Filtramos para solo dejar las comunas que tienen valor en hacinamiento
datos_prioridad_comunas_hacinamiento <- subset(dataframe_prioridad_comunas, dataframe_prioridad_comunas$Comuna %in% comunas_con_hacinamiento)
#Filtrar las de hacinamiento
datos_con_hacinamiento <- subset(regiones_y_comunas, regiones_y_comunas$comuna_nom %in% comunas_con_hacinamiento)
#Creanis un dataset con las comunas y su respectivo índice de hacinamiento y educación
indice_comunas_hacinamiento <- cbind(con_comunas$comuna_nom,con_comunas$hacinamiento)
indice_educacion_comuna <- cbind(con_comunas$comuna_nom,con_comunas$educ_sup_com)
carencia_servicios_bas <- cbind(con_comunas$comuna_nom,con_comunas$carente_ser_bas)Agregamos las columnas de hacinamiento y educación superior
#Re nombramos las columnas de los indices de hacinamiento y educación
colnames(indice_comunas_hacinamiento)[2] <- "Hacinamiento"
colnames(indice_comunas_hacinamiento)[1] <- "Comuna"
colnames(indice_educacion_comuna)[2] <- "Educacion_sup"
colnames(indice_educacion_comuna)[1] <- "Comuna"
colnames(carencia_servicios_bas)[2] <- "Servicios"
colnames(carencia_servicios_bas)[1] <- "Comuna"
#Agregamos las columnas con índice de hacinamiento y educación superior por comuna
datos_prioridad_comunas_hacinamiento <- merge(datos_prioridad_comunas_hacinamiento, indice_comunas_hacinamiento, by="Comuna", all.x=TRUE)
# Repetimos el proceso para agregar los datos de educación superior.
datos_prioridad_comunas_hacinamiento <- merge(datos_prioridad_comunas_hacinamiento, indice_educacion_comuna, by="Comuna", all.x=TRUE)
datos_prioridad_comunas_hacinamiento <- merge(datos_prioridad_comunas_hacinamiento, carencia_servicios_bas, by="Comuna", all.x=TRUE)
dactH=datos_prioridad_comunas_hacinamiento
dact=datos_prioridad_comunas_hacinamiento
dat=dias
Prio=dact$Prioridad
hacina = datos_prioridad_comunas_hacinamiento$Hacinamiento
hacina=datos_prioridad_comunas_hacinamiento$Hacinamiento
educa=datos_prioridad_comunas_hacinamiento$Educacion_sup
carencia_ser=educa=datos_prioridad_comunas_hacinamiento$Servicios
final_data <- datos_prioridad_comunas_hacinamientoSE ASIGNAN NÚMEROS A LAS PRIORIDADES
Se asignan números a las prioridades de las comunas haciendo agrupaciones
Prioaux=rep(0,length(Prio))
Prioaux[Prio=="Sin Prioridad"]=1
Prioaux[Prio=="Baja Prioridad"]=2
Prioaux[Prio=="Media Baja Prioridad"]=3
Prioaux[Prio=="Media Alta Prioridad"]=4
Prioaux[Prio=="Alta Prioridad"]=5
Prioaux=sort(Prioaux)
Prioaux[Prioaux==1]="Sin Prioridad"
Prioaux[Prioaux==2]="Baja Prioridad"
Prioaux[Prioaux==3]="Media Baja Prioridad"
Prioaux[Prioaux==4]="Media Alta Prioridad"
Prioaux[Prioaux==5]="Alta Prioridad"
datact=data.frame(
"prio"=Prioaux,
"comuna" = dact$Comuna
)SUAVIZADO
Se realiza el suavizado a los datos funcionales para una mejor obtención de gráficas
d.f=dact[,2:290]*1000
row.names(d.f)=dact$Comuna
fdact=fdata(d.f, argvals = dat$Dias, rangeval = c(1,1016))
#suavizando
out1<-optim.basis(fdact,type.CV=CV.S)
bsp<-create.bspline.basis(fdact$rangeval,nbasis=out1$numbasis.opt)
S.bsp <- S.basis(fdact$argvals, bsp)
sfda<-fdact
sfda$data<-fdact$data%*%S.bspGRAFICA VELOCIDAD DE CONTAGIO
x.d1 <- fdata.deriv(sfda)
plot(x.d1)
GRÁFICA CONTAGIOS
plot(abs(sfda))
#abline(h=0,lwd=2)
plot(t(dat),dact[2,2:290],type = "l")
for (i in 2:34){
lines(t(dat),dact[i,2:290],type = "l")}
split.fdata<-function(x,f,drop=FALSE,...){
if (!is.factor(f)) f<-factor(f)
nlev<-nlevels(f)
lev<-levels(f)
if (is.matrix(x$data)) x$data<-data.frame(x$data)
if (is.fdata(x)) {
out<-split(x$data,f,drop=drop,...)
}
for (i in 1:nlev) out[[lev[i]]]<-fdata(out[[lev[i]]],x$argvals,x$rangeval,x$names)
out
}TEST DE PERMUTACIÓN Se quiere ver si hay una diferencia significativa entre las comunas priorizadas y las no priorizadas
Prioaux[Prio=="Sin Prioridad"]=0
Prioaux[Prio=="Baja Prioridad"]=1
Prioaux[Prio=="Media Baja Prioridad"]=1
Prioaux[Prio=="Media Alta Prioridad"]=1
Prioaux[Prio=="Alta Prioridad"]=1
l1<-split.fdata(fdact,Prioaux)
plot(l1[[1]],col=1,ylim=c(0,10),lty="dotted")
lines(l1[[2]],col=2,lty="dotted")
lines(mean(l1[[1]]),lwd=2,col=1)
lines(mean(l1[[2]]),lwd=2,col=2)
dev.off()## null device
## 1temp.fd1=fdata2fd(l1[[1]])
temp.fd2=fdata2fd(l1[[2]])
#Test de permutación
dif.test=tperm.fd(temp.fd1,temp.fd2,q=0.1)
for (i in 1:length(dif.test$Tvals)){
if (dif.test$Tvals[i] >= dif.test$qval){
abline(v=dif.test$argvals[i])
}
}
dif.test=tperm.fd(temp.fd1,temp.fd2,q=0.1)
abline(v=50)
abline(v=88)
abline(v=326.7)
#abline(v=341.6)
#abline(v=380)
abline(v=422)
abline(v=599)
abline(v=615.5)
abline(v=713)
abline(v=805)GRÁFICA CONTAGIOS CON RANGOS DE DIFERENCIA SIGNIFICATIVA DE CONTAGIOS ENTRE COMUNAS PRIORIZADAS Y NO PRIORIZADAS
# Asignar colores según la prioridad
# Asignar colores según la prioridad
Prioaux[Prio=="Sin Prioridad"] = 0
Prioaux[Prio=="Baja Prioridad"] = 1
Prioaux[Prio=="Media Baja Prioridad"] = 1
Prioaux[Prio=="Media Alta Prioridad"] = 1
Prioaux[Prio=="Alta Prioridad"] = 1
colors <- ifelse(Prioaux == 0, "blue", "red")
# Plot con colores y leyenda
plot(abs(sfda), col = colors, main = "Contagions by communes with priority and without priority")
abline(v=50)
abline(v=88)
abline(v=326.7)
abline(v=422)
abline(v=599)
abline(v=615.5)
abline(v=713)
abline(v=805)
# Añadir línea punteada en la posición 280 y etiquetarla
abline(v=280, col="orange", lty=2)
text(280, max(abs(sfda)), "Vaccination", pos=4, col="green", cex=0.7)
# Añadir leyenda
legend("topright", legend=c("Non-prioritized", "Prioritized"), fill=c("blue", "red"))
GRÁFICA VELOCIDAD POR PRIORIDAD
Prioaux[Prio=="Sin Prioridad"]=0
Prioaux[Prio=="Baja Prioridad"]=1
Prioaux[Prio=="Media Baja Prioridad"]=1
Prioaux[Prio=="Media Alta Prioridad"]=1
Prioaux[Prio=="Alta Prioridad"]=1
l1<-split.fdata(x.d1,Prioaux)
#dev.off()
plot(l1[[1]],col=1,ylim=c(-0.4,0.4),lty="dotted")
lines(l1[[2]],col=2,lty="dotted")
lines(mean(l1[[1]]),lwd=2,col=1)
lines(mean(l1[[2]]),lwd=2,col=2)
TEST DE PERMUTACIÓN VELOCIDADES DE CONTAGIO
temp.fd1=fdata2fd(l1[[1]])
temp.fd2=fdata2fd(l1[[2]])
#Test de permutación
dif.test=tperm.fd(temp.fd1,temp.fd2,q=0.1)
for (i in 1:length(dif.test$Tvals)){
if (dif.test$Tvals[i] >= dif.test$qval){
abline(v=dif.test$argvals[i])
}
}
dif.test=tperm.fd(temp.fd1,temp.fd2,q=0.1)
abline(v=50);
abline(v=88);
abline(v=326.7);
#abline(v=341.6)
#abline(v=380)
abline(v=422);
abline(v=599);
abline(v=615.5);
abline(v=713);
abline(v=805);
Análisis con hacinamiento
# Asignar colores según la prioridad
Hacinamiento <- rep(0, length(hacina))
# hasta 0.01 - sin hacinamiento; ---> 0
Hacinamiento[hacina <= 0.01] <- 0
# de 0.01 a 0.04 - hacinamiento medio; ----> 1
Hacinamiento[hacina > 0.01] <- 1
# más de 0.04 - hacinamiento crítico. ----> 2
#Hacinamiento[hacina > 0.04] <- 2
colors <- ifelse(Hacinamiento == 0, "blue", "red")
# Plot con colores y leyenda
plot(abs(sfda), col = colors, main = "Contagios por comunas con diferencia de hacinamiento")
# Añadir leyenda
legend("topright",
legend = c("Sin Hacinamiento", "Hacinamiento Medio"), # Texto de la leyenda
fill = c("blue", "red")) # Colores correspondientes
# Añadir línea punteada en la posición 280 y etiquetarla
abline(v=280, col="orange", lty=2)
Velocidades de contagio por clasificación de hacinamiento
Hacinamiento <- rep(0, length(hacina))
# hasta 0.01 - sin hacinamiento; ---> 0
Hacinamiento[hacina <= 0.01] <- 0
# de 0.01 a 0.04 - hacinamiento medio; ----> 1
Hacinamiento[hacina > 0.01] <- 1
# más de 0.04 - hacinamiento crítico. ----> 2
#Hacinamiento[hacina > 0.04] <- 2
colors <- ifelse(Hacinamiento == 0, "blue", "red")
# Plot con colores y leyenda
plot(abs(x.d1), col = colors, main = "Contagios por comunas con diferencia de hacinamiento")
# Añadir leyenda
legend("topright",
legend = c("Sin Hacinamiento", "Hacinamiento Medio / critico"), # Texto de la leyenda
fill = c("blue", "red")) # Colores correspondientes
# Añadir línea punteada en la posición 280 y etiquetarla
abline(v=280, col="orange", lty=2)
Test de permutación hacinamiento y contagios
l1<-split.fdata(x.d1,Hacinamiento)
temp.fd1=fdata2fd(l1[[1]])
temp.fd2=fdata2fd(l1[[2]])
#Test de permutación
dif.test=tperm.fd(temp.fd1,temp.fd2,q=0.1)
for (i in 1:length(dif.test$Tvals)){
if (dif.test$Tvals[i] >= dif.test$qval){
abline(v=dif.test$argvals[i])
}
}
dif.test=tperm.fd(temp.fd1,temp.fd2,q=0.1) 
Diferencia de medias hacinamiento y velocidad de contagio
l1<-split.fdata(x.d1,Hacinamiento)
temp.fd1=fdata2fd(l1[[1]])
temp.fd2=fdata2fd(l1[[2]])
#Test de permutación
dif.test=tperm.fd(temp.fd1,temp.fd2,q=0.1)
for (i in 1:length(dif.test$Tvals)){
if (dif.test$Tvals[i] >= dif.test$qval){
abline(v=dif.test$argvals[i])
}
}
dif.test=tperm.fd(temp.fd1,temp.fd2,q=0.1) 
Intento prepración de datos
# Ejemplo simplificado de cómo podrías calcular el promedio de contagios para cada comuna
# Asumiendo que 'sfda' es una matriz con 27 filas (comunas) y 289 columnas (días)
promedio_contagios = apply(dact[2:200],1,mean)
datos_x.d1 <- x.d1$data
# Convierte los datos en una matriz
matriz_x.d1 <- as.matrix(datos_x.d1)
promedio_velocidad = apply(matriz_x.d1[,1:289],1,mean)
promedio_velocidad <- as.vector(promedio_velocidad)
#promedio_contagios <- apply(sfda, 1, mean) # Calcula el promedio para cada
carencia_ser = as.numeric(carencia_ser)
hacina = as.numeric(hacina)
modelo_regresion <- lm(promedio_velocidad ~ hacina + carencia_ser)
# Muestra el resumen del modelo
summary(modelo_regresion)##
## Call:
## lm(formula = promedio_velocidad ~ hacina + carencia_ser)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0012330 -0.0002249 0.0001166 0.0003217 0.0005625
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0001016 0.0001615 -0.629 0.535
## hacina 0.0103402 0.0046277 2.234 0.035 *
## carencia_ser 0.0017204 0.0037934 0.454 0.654
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0004325 on 24 degrees of freedom
## Multiple R-squared: 0.2061, Adjusted R-squared: 0.1399
## F-statistic: 3.115 on 2 and 24 DF, p-value: 0.0627Carencia servicios básicos
#promedio_contagios <- apply(sfda, 1, mean) # Calcula el promedio para cada
carencia_ser = as.numeric(carencia_ser)
modelo_regresion <- lm(promedio_contagios ~ carencia_ser)
# Muestra el resumen del modelo
summary(modelo_regresion)##
## Call:
## lm(formula = promedio_contagios ~ carencia_ser)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.513e-04 -6.558e-05 -5.857e-06 9.166e-05 2.316e-04
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.715e-04 3.757e-05 20.532 <2e-16 ***
## carencia_ser 1.500e-04 1.021e-03 0.147 0.884
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.000121 on 25 degrees of freedom
## Multiple R-squared: 0.0008634, Adjusted R-squared: -0.0391
## F-statistic: 0.0216 on 1 and 25 DF, p-value: 0.8843Regresión
Covarianza entre contagio y ducación superior, hacinamiento
# Extraemos las columnas relevantes de hacinamiento y contagios
indices_hacinamiento <- as.numeric(final_data$Hacinamiento)
indices_educacion <- as.numeric(final_data$Educacion_sup)
# Calculamos la covarianza entre los contagios y los índices de hacinamiento
cov_contagios_hacinamiento <- cov(sfda$data, indices_hacinamiento)
# Calculamos la covarianza entre los contagios y los índices de educación
cov_contagios_educacion <- cov(sfda$data, indices_educacion)
plot(cov_contagios_hacinamiento)
plot(cov_contagios_educacion)
Covarianza entre velocidad de contagio y ducación superior, hacinamiento
indices_hacinamiento <- as.numeric(final_data$Hacinamiento)
indices_educacion <- as.numeric(final_data$Educacion_sup)
# Calculamos la covarianza entre los contagios y los índices de hacinamiento
cov_contagios_hacinamiento <- cov(x.d1$data, indices_hacinamiento)
# Calculamos la covarianza entre los contagios y los índices de educación
cov_contagios_educacion <- cov(x.d1$data, indices_educacion)
# Imprimimos los resultados
plot(cov_contagios_hacinamiento)
plot(cov_contagios_educacion)
Covarianza entre la velocidad de contagio y los índices de hacinamiento
datos <- data.frame(Contagios = x.d1$data, Hacinamiento = indices_hacinamiento)
cov_mat <- cov(datos)
plot(cov_mat, main = "Matriz de Covarianza")
Primero necesitas transformar estos datos a un formato funcional si no lo están
contagios_fd <- smooth.basis(dat$Dias, t(as.matrix(final_data[,2:290])), fdParobj=fdPar(create.fourier.basis(range(dat$Dias), nbasis=20)))
# Extrayendo el vector de hacinamiento
hacinamiento <- as.numeric(hacina)
educacion <- as.numeric(educa)
# Configurando las covariables para la regresión (incluyendo un intercepto)
covariables <- list(intercepto = rep(1, nrow(final_data)), hacinamiento = hacinamiento)
conbasis <- create.constant.basis(range(dat$Dias))
betalist <- list(conbasis, create.fourier.basis(range(dat$Dias), nbasis=25))
# Ejecutando la regresión funcional
resultados <- fRegress.fd(contagios_fd$fd, covariables,betalist=betalist, wt=NULL,returnMatrix = FALSE,method=c('fRegress', 'model'),CVobs=1:N)
#std_regres <- fRegress.stderr(resultados,resultados$y2cMap,resultados$SigmaE,returnMatrix = FALSE,)
# Ejecutando la regresión funcional
# Extraer y visualizar resultados
beta_estimados <- resultados$betaest
summary(resultados)## Length Class Mode
## yfdobj 3 fd list
## xfdlist 2 -none- list
## betalist 2 -none- list
## betaestlist 2 -none- list
## yhatfdobj 3 fd list
## Cmat 676 -none- numeric
## Dmat 26 -none- numeric
## Cmatinv 676 -none- numeric
## wt 27 -none- numeric
## df 1 -none- logical
## y2cMap 0 -none- NULL
## SigmaE 0 -none- NULL
## betastderrlist 0 -none- NULL
## bvar 0 -none- NULL
## c2bMap 0 -none- NULLbeta_estimados_hacinamiento <- resultados$betaest
# Acceder a la base funcional y a los coeficientes para el hacinamiento
basis_hacina <- resultados$betaest[[2]]$fd$basis
coefs_hacina <- resultados$betaest[[2]]$fd$coefs
# Crear el objeto fd para el efecto del hacinamiento
fd_hacinamiento <- fd(coefs_hacina, basis_hacina)
summary(fd_hacinamiento)## Functional data object:
##
## Dimensions of the data:
## args
## reps
## funs
##
## Basis object:
##
## Type: fourier
##
## Range: 1 to 1016
##
## Number of basis functions: 25
##
## Period: 1015
##
## Coefficient matrix:## reps 1
## [1,] -8.211026e-02
## [2,] -1.581963e-01
## [3,] -3.371762e-02
## [4,] 4.168201e-02
## [5,] -5.700732e-02
## [6,] 2.815561e-01
## [7,] 2.185177e-01
## [8,] -1.010943e-01
## [9,] -1.842950e-01
## [10,] 1.728852e-01
## [11,] -1.025518e-01
## [12,] 5.884808e-02
## [13,] -1.116058e-03
## [14,] -2.387566e-01
## [15,] -1.230149e-01
## [16,] 8.453192e-02
## [17,] -2.467235e-01
## [18,] 1.563305e-02
## [19,] 1.203087e-01
## [20,] -2.245087e-01
## [21,] 7.460581e-03
## [22,] 8.608482e-05
## [23,] 4.023749e-05
## [24,] -3.889432e-05
## [25,] 9.070494e-05# Graficar el efecto del hacinamiento
plot(fd_hacinamiento, xlab="Dia", ylab="Efecto del indice de hacinamiento en los contagios")
## [1] "done"P_values de eefcto del hacinamiento en los contagios
coeficientes <- resultados$betaestlist[[2]]$fd$coefs
coeficientes <- rbind(resultados$betaestlist[[1]]$fd$coefs,coeficientes)
# Obtener el número de observaciones
N <- length(resultados$wt)
# Calcular la desviación estándar de los coeficientes
coef_std <- sqrt(diag(resultados$Cmatinv))
# Calcular el estadístico t para cada coeficiente
t_stats <- coeficientes / coef_std
# Calcular los p-valores para los coeficientes
p_values <- 2 * (1 - pt(abs(t_stats), df = N - length(coeficientes)))
# Mostrar los coeficientes, los t-stats y los p-valores
resultados_summary <- data.frame(coeficientes, coef_std, t_stats, p_values)
rownames(resultados_summary) <- names(coeficientes)
print(resultados_summary)## coeficientes coef_std t_stats p_values
## 1 9.486172e-04 0.01032241 9.189880e-02 0.9416593
## 2 -8.211026e-02 10.28588222 -7.982812e-03 0.9949181
## 3 -1.581963e-01 6.01929265 -2.628154e-02 0.9832725
## 4 -3.371762e-02 6.01929265 -5.601591e-03 0.9964340
## 5 4.168201e-02 6.01929265 6.924736e-03 0.9955916
## 6 -5.700732e-02 6.01929265 -9.470768e-03 0.9939709
## 7 2.815561e-01 6.01929265 4.677560e-02 0.9702434
## 8 2.185177e-01 6.01929265 3.630289e-02 0.9768990
## 9 -1.010943e-01 6.01929265 -1.679504e-02 0.9893089
## 10 -1.842950e-01 6.01929265 -3.061738e-02 0.9805145
## 11 1.728852e-01 6.01929265 2.872184e-02 0.9817201
## 12 -1.025518e-01 6.01929265 -1.703719e-02 0.9891548
## 13 5.884808e-02 6.01929265 9.776577e-03 0.9937762
## 14 -1.116058e-03 6.01929265 -1.854136e-04 0.9998820
## 15 -2.387566e-01 6.01929265 -3.966522e-02 0.9747616
## 16 -1.230149e-01 6.01929265 -2.043678e-02 0.9869914
## 17 8.453192e-02 6.01929680 1.404349e-02 0.9910602
## 18 -2.467235e-01 6.01928851 -4.098881e-02 0.9739203
## 19 1.563305e-02 6.01929265 2.597157e-03 0.9983466
## 20 1.203087e-01 6.01929265 1.998718e-02 0.9872775
## 21 -2.245087e-01 6.01929265 -3.729819e-02 0.9762662
## 22 7.460581e-03 6.01929265 1.239445e-03 0.9992109
## 23 8.608482e-05 6.01929265 1.430148e-05 0.9999909
## 24 4.023749e-05 6.01929265 6.684753e-06 0.9999957
## 25 -3.889432e-05 6.01929265 -6.461610e-06 0.9999959
## 26 9.070494e-05 6.01929265 1.506904e-05 0.9999904beta_estimados_hacinamiento <- resultados$betaest
# Extraer la matriz de penalización y varianzas-covarianzas
Cmatinv <- solve(resultados$Cmat)
sigmaHat <- diag(Cmatinv) / length(dat$Dias)
# Obtener los errores estándar de los coeficientes
stderrlist <- lapply(resultados$betaest, function(beta) {
basisobj <- beta$fd$basis
nbasis <- basisobj$nbasis
sqrt(diag(Cmatinv[1:nbasis, 1:nbasis]))
})
# Calcular los p-values
p_values_hacina <- sapply(2:length(stderrlist), function(i) {
coefs <- resultados$betaest[[i]]$fd$coefs
stderr <- stderrlist[[i]]
2 * (1 - pnorm(abs(coefs / stderr)))
})
# Mostrar los p-values
#print(p_values)
coefs_hacina <- resultados$betaestlist[[2]]$fd$coefs
resultados_summary_hacina <- data.frame(coefs_hacina, p_values_hacina)
names = c("Coeficiente", "p_values")
rownames(resultados_summary_hacina) <- names(names)
print(resultados_summary_hacina)## coefs_hacina p_values_hacina
## 1 -8.211026e-02 1.776357e-15
## 2 -1.581963e-01 9.877291e-01
## 3 -3.371762e-02 9.955306e-01
## 4 4.168201e-02 9.944749e-01
## 5 -5.700732e-02 9.924435e-01
## 6 2.815561e-01 9.626921e-01
## 7 2.185177e-01 9.710408e-01
## 8 -1.010943e-01 9.866001e-01
## 9 -1.842950e-01 9.755747e-01
## 10 1.728852e-01 9.770864e-01
## 11 -1.025518e-01 9.864069e-01
## 12 5.884808e-02 9.921995e-01
## 13 -1.116058e-03 9.998521e-01
## 14 -2.387566e-01 9.683600e-01
## 15 -1.230149e-01 9.836949e-01
## 16 8.453192e-02 9.887953e-01
## 17 -2.467235e-01 9.673049e-01
## 18 1.563305e-02 9.979278e-01
## 19 1.203087e-01 9.840536e-01
## 20 -2.245087e-01 9.702473e-01
## 21 7.460581e-03 9.990111e-01
## 22 8.608482e-05 9.999886e-01
## 23 4.023749e-05 9.999947e-01
## 24 -3.889432e-05 9.999948e-01
## 25 9.070494e-05 9.999880e-01Interpretación de gráfico:
Eje X (Día): Representa el tiempo. Va de 0 a 1016.
Eje Y (Efecto del Hacinamiento en Contagios): Muestra la magnitud y la dirección del efecto que el hacinamiento tiene en los contagios. Los valores positivos indican un aumento en los contagios con el aumento del hacinamiento, y los valores negativos indican una disminución en los contagios con el aumento del hacinamiento.
La línea en el gráfico: Es la función estimada que describe cómo cambia el efecto del hacinamiento en los contagios a lo largo del tiempo. La forma de la línea es una combinación de funciones sinusoidales y cosinusoidales debido a que se utilizó una base de Fourier para estimar estos efectos.
Variación temporal: Hay una variación clara a lo largo del tiempo en el efecto del hacinamiento sobre los contagios, lo que indica que este efecto no es constante.
Pico y valles: Parece haber un pico prominente que indica que en cierto punto, el efecto del hacinamiento sobre los contagios fue positivo y bastante significativo. Esto es seguido por un valle donde el efecto se vuelve negativo.
Cruce por cero: El efecto cruza la línea de cero varias veces, lo que indica que hay momentos en los cuales el hacinamiento no tiene un efecto estadísticamente significativo sobre los contagios.
Regresión velocidad de contagios y hacinamiento
velocidad_fd <- smooth.basis(dat$Dias, t(as.matrix(x.d1$data[,1:289])), fdParobj=fdPar(create.fourier.basis(range(dat$Dias), nbasis=20)))
# Extrayendo el vector de hacinamiento
hacinamiento <- as.numeric(hacina)
educacion <- as.numeric(educa)
# Configurando las covariables para la regresión (incluyendo un intercepto)
covariables_velocidad <- list(intercepto = rep(1, nrow(final_data)), hacinamiento = hacinamiento )
conbasis <- create.constant.basis(range(dat$Dias))
betalist <- list(conbasis, create.fourier.basis(range(dat$Dias), nbasis=20))
# Ejecutando la regresión funcional
regresion_velocidad_hacinamiento <- fRegress(velocidad_fd$fd, covariables_velocidad, betalist=betalist)
# Ejecutando la regresión funcional
# Extraer y visualizar resultados
beta_estimados_velocidad_hacinamiento <- regresion_velocidad_hacinamiento$betaest
summary(regresion_velocidad_hacinamiento)## Length Class Mode
## yfdobj 3 fd list
## xfdlist 2 -none- list
## betalist 2 -none- list
## betaestlist 2 -none- list
## yhatfdobj 3 fd list
## Cmat 484 -none- numeric
## Dmat 22 -none- numeric
## Cmatinv 484 -none- numeric
## wt 27 -none- numeric
## df 1 -none- logical
## y2cMap 0 -none- NULL
## SigmaE 0 -none- NULL
## betastderrlist 0 -none- NULL
## bvar 0 -none- NULL
## c2bMap 0 -none- NULL# Acceder a la base funcional y a los coeficientes para el hacinamiento
basis_reg_velo_hacina <- regresion_velocidad_hacinamiento$betaest[[2]]$fd$basis
coefs_velo_hacina <- regresion_velocidad_hacinamiento$betaest[[2]]$fd$coefs
# Crear el objeto fd para el efecto del hacinamiento
fd_velo_hacinamiento <- fd(coefs_velo_hacina, basis_reg_velo_hacina)
summary(regresion_velocidad_hacinamiento)## Length Class Mode
## yfdobj 3 fd list
## xfdlist 2 -none- list
## betalist 2 -none- list
## betaestlist 2 -none- list
## yhatfdobj 3 fd list
## Cmat 484 -none- numeric
## Dmat 22 -none- numeric
## Cmatinv 484 -none- numeric
## wt 27 -none- numeric
## df 1 -none- logical
## y2cMap 0 -none- NULL
## SigmaE 0 -none- NULL
## betastderrlist 0 -none- NULL
## bvar 0 -none- NULL
## c2bMap 0 -none- NULL# Graficar el efecto del hacinamiento
plot(fd_velo_hacinamiento, xlab="Dia", ylab="Efecto del indice de hacinamiento en la velocidad de contagio")
## [1] "done"Tabla efecto del hacinamiento en la velocidad de contagio
coeficientes_velo_hacina <- regresion_velocidad_hacinamiento$betaestlist[[2]]$fd$coefs
coeficientes_velo_hacina <- rbind(regresion_velocidad_hacinamiento$betaestlist[[1]]$fd$coefs,coeficientes_velo_hacina)
# Obtener el número de observaciones
N <- length(regresion_velocidad_hacinamiento$wt)
# Calcular la desviación estándar de los coeficientes
coef_std_velo_hacina <- sqrt(diag(regresion_velocidad_hacinamiento$Cmatinv))
# Calcular el estadístico t para cada coeficiente
t_stats_velo_hacina <- coeficientes_velo_hacina / coef_std_velo_hacina
# Calcular los p-valores para los coeficientes
# Obtener el número de coeficientes
num_coeficientes <- ncol(coeficientes_velo_hacina)
# Calcular los p-valores para los coeficientes con los grados de libertad correctos
p_values_velo_hacina <- 2 * (1 - pt(abs(t_stats_velo_hacina), df = N - num_coeficientes))
# Mostrar los coeficientes, los t-stats y los p-valores
resultados_summary_velo_hacina <- data.frame(coeficientes_velo_hacina, coef_std_velo_hacina, t_stats_velo_hacina, p_values_velo_hacina)
rownames(resultados_summary_velo_hacina) <- names(coeficientes_velo_hacina)
print(resultados_summary_velo_hacina)## coeficientes_velo_hacina coef_std_velo_hacina t_stats_velo_hacina
## 1 -8.771943e-05 0.01032241 -0.008497959
## 2 2.189059e-01 10.28588222 0.021282167
## 3 4.185541e-01 6.01929265 0.069535428
## 4 -7.710826e-01 6.01929265 -0.128101864
## 5 1.038467e+00 6.01929265 0.172523049
## 6 7.177524e-01 6.01929265 0.119241978
## 7 -3.740021e+00 6.01929265 -0.621338951
## 8 5.411947e+00 6.01929265 0.899100142
## 9 4.736249e+00 6.01929265 0.786844845
## 10 -2.361175e+00 6.01929265 -0.392267794
## 11 3.130939e+00 6.01929265 0.520150569
## 12 5.430371e+00 6.01929265 0.902160990
## 13 -2.092536e-01 6.01929265 -0.034763825
## 14 2.211494e+00 6.01929265 0.367400986
## 15 4.951017e+00 6.01929265 0.822524694
## 16 -1.035323e+01 6.01929265 -1.720008172
## 17 1.184726e+01 6.01929680 1.968213948
## 18 4.192892e+00 6.01928851 0.696576025
## 19 -6.967792e+00 6.01929265 -1.157576577
## 20 9.487723e-01 6.01929265 0.157621900
## 21 -5.859088e-01 6.01929265 -0.097338474
## 22 -1.375190e+01 6.01929265 -2.284637139
## p_values_velo_hacina
## 1 0.99328456
## 2 0.98318306
## 3 0.94509549
## 4 0.89905408
## 5 0.86436230
## 6 0.90600022
## 7 0.53978657
## 8 0.37685116
## 9 0.43848649
## 10 0.69805758
## 11 0.60735972
## 12 0.37525347
## 13 0.97253359
## 14 0.71629367
## 15 0.41825838
## 16 0.09731452
## 17 0.05979308
## 18 0.49224937
## 19 0.25756084
## 20 0.87597250
## 21 0.92320409
## 22 0.03073908a
library(fda)
# Configuración inicial
velocidad_fd <- smooth.basis(dat$Dias, t(as.matrix(x.d1$data[,1:289])), fdParobj=fdPar(create.fourier.basis(range(dat$Dias), nbasis=20)))
hacinamiento <- as.numeric(hacina)
covariables_velocidad <- list(intercepto = rep(1, nrow(final_data)), hacinamiento = hacinamiento)
conbasis <- create.constant.basis(range(dat$Dias))
betalist <- list(conbasis, create.fourier.basis(range(dat$Dias), nbasis=20))
# Ejecutar la regresión funcional
regresion_velocidad_hacinamiento <- fRegress(velocidad_fd$fd, covariables_velocidad, betalist=betalist)
# Extraer y visualizar resultados
beta_estimados_velocidad_hacinamiento <- regresion_velocidad_hacinamiento$betaest
# Extraer la matriz de penalización y varianzas-covarianzas
Cmatinv <- solve(regresion_velocidad_hacinamiento$Cmat)
sigmaHat <- diag(Cmatinv) / length(dat$Dias)
# Obtener los errores estándar de los coeficientes
stderrlist <- lapply(regresion_velocidad_hacinamiento$betaest, function(beta) {
basisobj <- beta$fd$basis
nbasis <- basisobj$nbasis
sqrt(diag(Cmatinv[1:nbasis, 1:nbasis]))
})
# Calcular los p-values
p_values <- sapply(2:length(stderrlist), function(i) {
coefs <- regresion_velocidad_hacinamiento$betaest[[i]]$fd$coefs
stderr <- stderrlist[[i]]
2 * (1 - pnorm(abs(coefs / stderr)))
})
# Mostrar los p-values
print(p_values)## [,1]
## [1,] 0.00000000
## [2,] 0.96754136
## [3,] 0.89806836
## [4,] 0.86302634
## [5,] 0.90508365
## [6,] 0.53437663
## [7,] 0.36859932
## [8,] 0.43137270
## [9,] 0.69486036
## [10,] 0.60295863
## [11,] 0.36697135
## [12,] 0.97226807
## [13,] 0.71331994
## [14,] 0.41077834
## [15,] 0.08543096
## [16,] 0.04904327
## [17,] 0.48606876
## [18,] 0.24703655
## [19,] 0.87475475
## [20,] 0.92245760
## [21,] 0.02233411coefs_velo_hacina <- regresion_velocidad_hacinamiento$betaestlist[[2]]$fd$coefs
resultados_summary_velo_hacina <- data.frame(coefs_velo_hacina, p_values)
names = c("Coeficiente", "p_values")
rownames(resultados_summary_velo_hacina) <- names(names)
print(resultados_summary_velo_hacina)## coefs_velo_hacina p_values
## 1 0.2189059 0.00000000
## 2 0.4185541 0.96754136
## 3 -0.7710826 0.89806836
## 4 1.0384667 0.86302634
## 5 0.7177524 0.90508365
## 6 -3.7400210 0.53437663
## 7 5.4119469 0.36859932
## 8 4.7362494 0.43137270
## 9 -2.3611747 0.69486036
## 10 3.1309385 0.60295863
## 11 5.4303710 0.36697135
## 12 -0.2092536 0.97226807
## 13 2.2114941 0.71331994
## 14 4.9510168 0.41077834
## 15 -10.3532326 0.08543096
## 16 11.8472639 0.04904327
## 17 4.1928921 0.48606876
## 18 -6.9677922 0.24703655
## 19 0.9487723 0.87475475
## 20 -0.5859088 0.92245760
## 21 -13.7518996 0.02233411Regresión contagios y nivel de educación
contagios_fd <- smooth.basis(dat$Dias, t(as.matrix(final_data[,2:290])), fdParobj=fdPar(create.fourier.basis(range(dat$Dias), nbasis=20)))
# Extrayendo el vector de hacinamiento
educacion <- as.numeric(educa)
# Configurando las covariables para la regresión (incluyendo un intercepto)
covariables2 <- list(intercepto = rep(1, nrow(final_data)), educacion = educacion )
conbasis <- create.constant.basis(range(dat$Dias))
betalist <- list(conbasis, create.fourier.basis(range(dat$Dias), nbasis=20))
# Ejecutando la regresión funcional
resultados_educ <- fRegress(contagios_fd$fd, covariables2, betalist=betalist)
# Ejecutando la regresión funcional
# Extraer y visualizar resultados
beta_estimados_educ <- resultados$betaest
summary(resultados_educ)## Length Class Mode
## yfdobj 3 fd list
## xfdlist 2 -none- list
## betalist 2 -none- list
## betaestlist 2 -none- list
## yhatfdobj 3 fd list
## Cmat 484 -none- numeric
## Dmat 22 -none- numeric
## Cmatinv 484 -none- numeric
## wt 27 -none- numeric
## df 1 -none- logical
## y2cMap 0 -none- NULL
## SigmaE 0 -none- NULL
## betastderrlist 0 -none- NULL
## bvar 0 -none- NULL
## c2bMap 0 -none- NULL# Acceder a la base funcional y a los coeficientes para el hacinamiento
basis_educ <- resultados_educ$betaest[[2]]$fd$basis
coefs_educ <- resultados_educ$betaest[[2]]$fd$coefs
# Crear el objeto fd para el efecto del hacinamiento
fd_hacinamiento <- fd(coefs_educ, basis_educ)
# Graficar el efecto del hacinamiento
plot(fd_hacinamiento, xlab="Dia", ylab="Efecto del nivel de educacion en los contagios")
## [1] "done"P_values para nivel de educación
contagios_fd <- smooth.basis(dat$Dias, t(as.matrix(final_data[,2:290])), fdParobj=fdPar(create.fourier.basis(range(dat$Dias), nbasis=20)))
coeficientes_educ <- resultados_educ$betaestlist[[2]]$fd$coefs
coeficientes_educ2 <- rbind(resultados_educ$betaestlist[[1]]$fd$coefs,coeficientes_educ)
# Obtener el número de observaciones
N <- length(resultados_educ$wt)
# Calcular la desviación estándar de los coeficientes
coef_std_educ <- sqrt(diag(resultados_educ$Cmatinv))
# Calcular el estadístico t para cada coeficiente
t_stats_educ <- coeficientes_educ2 / coef_std_educ
# Calcular los p-valores para los coeficientes
p_values_educ <- 2 * (1 - pt(abs(t_stats_educ), df = N - length(coeficientes_educ2)))
# Mostrar los coeficientes, los t-stats y los p-valores
resultados_educ_summary <- data.frame(coeficientes_educ2, coef_std_educ, t_stats_educ, p_values_educ)
rownames(resultados_educ_summary) <- names(coeficientes_educ2)
print(resultados_educ_summary)## coeficientes_educ2 coef_std_educ t_stats_educ p_values_educ
## 1 0.0008958571 0.009743838 0.091940894 0.9303150
## 2 -0.0155041616 8.431491410 -0.001838840 0.9986039
## 3 -0.1779934720 5.227083735 -0.034052156 0.9741531
## 4 -0.0235982267 5.227083735 -0.004514607 0.9965725
## 5 0.0281403387 5.227083735 0.005383564 0.9959128
## 6 -0.0729308110 5.227083735 -0.013952486 0.9894075
## 7 0.2277096515 5.227083735 0.043563421 0.9669386
## 8 0.1781161929 5.227083735 0.034075634 0.9741353
## 9 -0.0911380464 5.227083735 -0.017435735 0.9867634
## 10 -0.1442532617 5.227083735 -0.027597274 0.9790510
## 11 0.1401416193 5.227083735 0.026810670 0.9796479
## 12 -0.0866762664 5.227083735 -0.016582146 0.9874113
## 13 0.0376793825 5.227083735 0.007208490 0.9945273
## 14 0.0149493877 5.227083735 0.002859986 0.9978287
## 15 -0.2053881709 5.227083735 -0.039293071 0.9701774
## 16 -0.1042886796 5.227083735 -0.019951599 0.9848537
## 17 0.0754336063 5.227087332 0.014431289 0.9890440
## 18 -0.2111672391 5.227080137 -0.040398699 0.9693388
## 19 0.0215730571 5.227083735 0.004127169 0.9968666
## 20 0.1122382568 5.227083735 0.021472443 0.9836993
## 21 -0.1868058769 5.227083735 -0.035738069 0.9728741
## 22 0.0074251156 5.227083735 0.001420508 0.9989215Regresión velocidad de contagio y educación
velocidad_fd <- smooth.basis(dat$Dias, t(as.matrix(x.d1$data[,1:289])), fdParobj=fdPar(create.fourier.basis(range(dat$Dias), nbasis=20)))
# Extrayendo el vector de hacinamiento
hacinamiento <- as.numeric(hacina)
educacion <- as.numeric(educa)
# Configurando las covariables para la regresión (incluyendo un intercepto)
covariables_velocidad_educ <- list(intercepto = rep(1, nrow(final_data)), educacion = educacion )
conbasis <- create.constant.basis(range(dat$Dias))
betalist <- list(conbasis, create.fourier.basis(range(dat$Dias), nbasis=20))
# Ejecutando la regresión funcional
regresion_velocidad_educ <- fRegress(velocidad_fd$fd, covariables_velocidad_educ, betalist=betalist)
# Ejecutando la regresión funcional
# Extraer y visualizar resultados
beta_estimados_velocidad_educ <- regresion_velocidad_educ$betaest
summary(regresion_velocidad_educ)## Length Class Mode
## yfdobj 3 fd list
## xfdlist 2 -none- list
## betalist 2 -none- list
## betaestlist 2 -none- list
## yhatfdobj 3 fd list
## Cmat 484 -none- numeric
## Dmat 22 -none- numeric
## Cmatinv 484 -none- numeric
## wt 27 -none- numeric
## df 1 -none- logical
## y2cMap 0 -none- NULL
## SigmaE 0 -none- NULL
## betastderrlist 0 -none- NULL
## bvar 0 -none- NULL
## c2bMap 0 -none- NULL# Acceder a la base funcional y a los coeficientes para el hacinamiento
basis_reg_velo_educ <- regresion_velocidad_educ$betaest[[2]]$fd$basis
coefs_velo_educ <- regresion_velocidad_educ$betaest[[2]]$fd$coefs
# Crear el objeto fd para el efecto del hacinamiento
fd_velo_educ <- fd(coefs_velo_educ, basis_reg_velo_educ)
summary(regresion_velocidad_educ)## Length Class Mode
## yfdobj 3 fd list
## xfdlist 2 -none- list
## betalist 2 -none- list
## betaestlist 2 -none- list
## yhatfdobj 3 fd list
## Cmat 484 -none- numeric
## Dmat 22 -none- numeric
## Cmatinv 484 -none- numeric
## wt 27 -none- numeric
## df 1 -none- logical
## y2cMap 0 -none- NULL
## SigmaE 0 -none- NULL
## betastderrlist 0 -none- NULL
## bvar 0 -none- NULL
## c2bMap 0 -none- NULL# Graficar el efecto del hacinamiento
plot(fd_velo_educ, xlab="Dia", ylab="Efecto del nivel de educacion en la velocidad de contagio")
## [1] "done"P_values relación entre velocidad de contagio y nivel de educación
# Extraer y visualizar resultados
beta_estimados_velocidad_educ <- regresion_velocidad_educ$betaest
# Extraer la matriz de penalización y varianzas-covarianzas
Cmatinv <- solve(regresion_velocidad_educ$Cmat)
sigmaHat <- diag(Cmatinv) / length(dat$Dias)
# Obtener los errores estándar de los coeficientes
stderrlist <- lapply(regresion_velocidad_educ$betaest, function(beta) {
basisobj <- beta$fd$basis
nbasis <- basisobj$nbasis
sqrt(diag(Cmatinv[1:nbasis, 1:nbasis]))
})
# Calcular los p-values
p_values_educ <- sapply(2:length(stderrlist), function(i) {
coefs <- regresion_velocidad_educ$betaest[[i]]$fd$coefs
stderr <- stderrlist[[i]]
2 * (1 - pnorm(abs(coefs / stderr)))
})
# Mostrar los p-values
#print(p_values_educ)
coefs_velo_educ <- regresion_velocidad_educ$betaest[[2]]$fd$coefs
resultados_summary_velo_educ <- data.frame(coefs_velo_educ, p_values_educ)
names = c("Coeficiente", "p_values")
rownames(resultados_summary_velo_educ) <- names(names)
print(resultados_summary_velo_educ)## coefs_velo_educ p_values_educ
## 1 0.1123137 0.00000000
## 2 0.3179861 0.96991563
## 3 -0.9513312 0.85558244
## 4 1.1739810 0.82229374
## 5 0.5029171 0.92335085
## 6 -3.0519515 0.55930596
## 7 4.3793350 0.40213398
## 8 3.7051000 0.47843161
## 9 -2.1318438 0.68338693
## 10 2.6367954 0.61394613
## 11 4.4124841 0.39858117
## 12 -0.7705106 0.88281049
## 13 1.4248820 0.78516380
## 14 4.2019234 0.42146901
## 15 -8.9117308 0.08821037
## 16 10.1537854 0.05207255
## 17 3.7260894 0.47594330
## 18 -6.4643109 0.21619983
## 19 1.2461219 0.81157307
## 20 -0.5537728 0.91562760
## 21 -11.4631111 0.02830577Analizando datos antes de la vacunación:
contagios_fd_prevac <- smooth.basis(dat$Dias[1:90], t(as.matrix(x.d1$data[,2:91])), fdParobj=fdPar(create.fourier.basis(range(dat$Dias[1:90]), nbasis=20)))
# Configurando las covariables para la regresión (incluyendo un intercepto)
covariables2 <- list(intercepto = rep(1, nrow(final_data[,2:91])), hacinamiento = hacinamiento )
conbasis <- create.constant.basis(range(dat$Dias[1:90]))
betalist_prevac <- list(conbasis, create.fourier.basis(range(dat$Dias[1:90]), nbasis=20))
# Ejecutando la regresión funcional
resultados_pre_vac <- fRegress(contagios_fd_prevac$fd, covariables2, betalist=betalist_prevac)
# Ejecutando la regresión funcional
# Extraer y visualizar resultados
beta_estimados_educ_prevac <- resultados_pre_vac$betaest
summary(resultados_pre_vac)## Length Class Mode
## yfdobj 3 fd list
## xfdlist 2 -none- list
## betalist 2 -none- list
## betaestlist 2 -none- list
## yhatfdobj 3 fd list
## Cmat 484 -none- numeric
## Dmat 22 -none- numeric
## Cmatinv 484 -none- numeric
## wt 27 -none- numeric
## df 1 -none- logical
## y2cMap 0 -none- NULL
## SigmaE 0 -none- NULL
## betastderrlist 0 -none- NULL
## bvar 0 -none- NULL
## c2bMap 0 -none- NULL# Acceder a la base funcional y a los coeficientes para el hacinamiento
basis_educ_pre_vac <- resultados_pre_vac$betaest[[2]]$fd$basis
coefs_educ_pre_vac <- resultados_pre_vac$betaest[[2]]$fd$coefs
# Crear el objeto fd para el efecto del hacinamiento
fd_educa_pre_vac <- fd(coefs_educ_pre_vac, basis_educ_pre_vac)
# Graficar el efecto del hacinamiento
plot(fd_educa_pre_vac, xlab="Dia", ylab="Efecto del nivel de educacion en los contagios")
## [1] "done"P_values pre-vacunación
# Extraer la matriz de penalización y varianzas-covarianzas
Cmatinv <- solve(resultados_pre_vac$Cmat)
sigmaHat <- diag(Cmatinv) / length(dat$Dias)
# Obtener los errores estándar de los coeficientes
stderrlist <- lapply(resultados_pre_vac$betaest, function(beta) {
basisobj <- beta$fd$basis
nbasis <- basisobj$nbasis
sqrt(diag(Cmatinv[1:nbasis, 1:nbasis]))
})
# Calcular los p-values
p_values_prevac <- sapply(2:length(stderrlist), function(i) {
coefs <- resultados_pre_vac$betaest[[i]]$fd$coefs
stderr <- stderrlist[[i]]
2 * (1 - pnorm(abs(coefs / stderr)))
})
# Mostrar los p-values
#print(p_values_educ)
coefs_velo_prevac <- resultados_pre_vac$betaest[[2]]$fd$coefs
resultados_summary_velo_prevac <- data.frame(coefs_velo_prevac, p_values_prevac)
names = c("Coeficiente", "p_values")
rownames(resultados_summary_velo_prevac) <- names(names)
print(resultados_summary_velo_prevac)## coefs_velo_prevac p_values_prevac
## 1 0.18113252 0.0000000
## 2 -0.37546376 0.9708814
## 3 5.46727096 0.3637244
## 4 7.03373231 0.2425925
## 5 1.06912139 0.8590247
## 6 1.86667725 0.7564730
## 7 -8.01008029 0.1832764
## 8 -7.00920942 0.2442387
## 9 -2.05235217 0.7331319
## 10 -0.32609251 0.9567961
## 11 5.23933823 0.3840686
## 12 4.36474463 0.4683748
## 13 -1.59498426 0.7910259
## 14 -1.98455119 0.7416279
## 15 -3.98934688 0.5074842
## 16 -1.13798367 0.8500489
## 17 3.09752437 0.6068332
## 18 1.71018252 0.7763207
## 19 -1.03657662 0.8632732
## 20 -1.26615923 0.8333944
## 21 -0.03600609 0.9952273Relación entre hacinamiento y educación
modelo <- lm(educa ~ hacina, data = dat)
# Mostrar un resumen del modelo
summary(modelo)##
## Call:
## lm(formula = educa ~ hacina, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.026897 -0.013621 -0.006897 0.008103 0.079741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.020172 0.007499 2.690 0.0125 *
## hacina 0.336207 0.234539 1.433 0.1641
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0228 on 25 degrees of freedom
## Multiple R-squared: 0.07595, Adjusted R-squared: 0.03899
## F-statistic: 2.055 on 1 and 25 DF, p-value: 0.1641educa = as.numeric(educa)
modelo2 <- lm(hacina ~ educa, data = dat)
# Mostrar un resumen del modelo
summary(modelo2)##
## Call:
## lm(formula = hacina ~ educa, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.030695 -0.014084 -0.000695 0.008341 0.043823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.019400 0.005802 3.343 0.00261 **
## educa 0.225908 0.157594 1.433 0.16411
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01869 on 25 degrees of freedom
## Multiple R-squared: 0.07595, Adjusted R-squared: 0.03899
## F-statistic: 2.055 on 1 and 25 DF, p-value: 0.1641