HPAI-risk-factors
Alfredo Acosta
2025-03-06
1 Libraries
library(tidyverse); library(lubridate); library(epiDisplay)## Warning: package 'forcats' was built under R version 4.3.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors## Warning: package 'epiDisplay' was built under R version 4.3.3## Loading required package: foreign
## Loading required package: survival
## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
##
## Loading required package: nnet
##
## Attaching package: 'epiDisplay'
##
## The following object is masked from 'package:ggplot2':
##
## alphalibrary(dplyr);library(sjmisc); library(ResourceSelection)##
## Attaching package: 'sjmisc'
##
## The following object is masked from 'package:purrr':
##
## is_empty
##
## The following object is masked from 'package:tidyr':
##
## replace_na
##
## The following object is masked from 'package:tibble':
##
## add_case## Warning: package 'ResourceSelection' was built under R version 4.3.3## ResourceSelection 0.3-6 2023-06-27library(modEvA); library(lubridate)## Warning: package 'modEvA' was built under R version 4.3.3#WD
2 Quick analysis
3 factors
v$fcaso <- factor(v$caso_o_notifi, levels = c(0,1), labels = c("Control","Caso"))
table(v$fcaso)##
## Control Caso
## 124 334 Correlation between numeric
library(PerformanceAnalytics)## Warning: package 'PerformanceAnalytics' was built under R version 4.3.3## Loading required package: xts## Warning: package 'xts' was built under R version 4.3.3## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## #################################################################################
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, last##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legendchart.Correlation(v[,c(13,16,25)])
colnames(v)## [1] "orden" "caso_o_notifi"
## [3] "provincia" "cantón"
## [5] "parroquia" "predio"
## [7] "t_explotación" "Nuevo_bioseguridad"
## [9] "Nuevo_tipo" "notificador"
## [11] "f_1_ave_enferma" "f_notificación"
## [13] "Nuevo_dias_entre_enf_notifi" "Nuevo_fecha 1raenf_notifi"
## [15] "f_1_visita" "Nuevo_dias_entre_notifi_visita"
## [17] "síndrome_presuntivo" "patología"
## [19] "especie_final" "f_elaboración"
## [21] "f_ingreso_sizse" "f_cierre"
## [23] "Nuevo_fecha_enf_cierre" "det_diagnóstico"
## [25] "poblacion" "Log_poblacion"
## [27] "Nuevo_tipo_comercio" "ano"
## [29] "definitivo2" "fcaso"str(v)## tibble [157 × 30] (S3: tbl_df/tbl/data.frame)
## $ orden : num [1:157] 116 117 118 119 120 121 122 123 124 125 ...
## $ caso_o_notifi : num [1:157] 0 0 0 0 0 0 0 0 0 1 ...
## $ provincia : chr [1:157] "LOJA" "ZAMORA CHINCHIPE" "LOJA" "LOJA" ...
## $ cantón : chr [1:157] "GONZANAMA" "ZAMORA" "OLMEDO" "LOJA" ...
## $ parroquia : chr [1:157] "GONZANAMA" "CUMBARATZA" "LA TINGUE" "VALLE" ...
## $ predio : chr [1:157] "EL GUATO" "San Agustin" "Tambara" "CC de Fauna Silvestre Orillas del Zamora" ...
## $ t_explotación : chr [1:157] "Incubación" "Aves de Traspatio" "Aves de Traspatio" "Aves Silvestres" ...
## $ Nuevo_bioseguridad : chr [1:157] "no" "no" "no" "si" ...
## $ Nuevo_tipo : chr [1:157] "TRASPATIO" "TRASPATIO" "TRASPATIO" "SILVESTRES" ...
## $ notificador : chr [1:157] "PROPIETARIO" "PROPIETARIO" "PROPIETARIO" "PROPIETARIO" ...
## $ f_1_ave_enferma : POSIXct[1:157], format: "2022-01-15" "2022-04-28" ...
## $ f_notificación : POSIXct[1:157], format: "2022-03-07" "2022-05-09" ...
## $ Nuevo_dias_entre_enf_notifi : num [1:157] 51 11 4 16 7 6 5 5 102 11 ...
## $ Nuevo_fecha 1raenf_notifi : chr [1:157] "> 1 semana" "> 1 semana" "1 semana" "> 1 semana" ...
## $ f_1_visita : POSIXct[1:157], format: "2022-03-07" "2022-05-09" ...
## $ Nuevo_dias_entre_notifi_visita: num [1:157] 0 0 1 0 0 0 1 0 0 0 ...
## $ síndrome_presuntivo : chr [1:157] "Respiratorio Agudo" "Respiratorio Crónico" "Respiratorio Crónico" "Gastroentérico" ...
## $ patología : chr [1:157] "ENFERMEDADES RESPIRATORIAS" "Enfermedad Respiratoria" "Respiratoria" "Digestiva" ...
## $ especie_final : chr [1:157] "NA" "PATOS" "NA" "AVES SILVESTRES" ...
## $ f_elaboración : POSIXct[1:157], format: "2022-03-08" "2022-05-09" ...
## $ f_ingreso_sizse : POSIXct[1:157], format: "2022-03-08 08:48:21" "2022-05-09 16:11:00" ...
## $ f_cierre : POSIXct[1:157], format: "2022-04-12" "2022-06-07" ...
## $ Nuevo_fecha_enf_cierre : num [1:157] 87 40 25 44 30 56 39 38 205 219 ...
## $ det_diagnóstico : chr [1:157] "Infección por Virus de Influenza Aviar" "Infección por Virus de Influenza Aviar" "Infección por Virus de Influenza Aviar" "Infección por Virus de Influenza Aviar" ...
## $ poblacion : num [1:157] 70 33 40 1 15 ...
## $ Log_poblacion : chr [1:157] "2" "2" "2" "0" ...
## $ Nuevo_tipo_comercio : chr [1:157] "TRASPATIO <500" "TRASPATIO <500" "TRASPATIO <500" "SILVESTRE" ...
## $ ano : num [1:157] 2022 2022 2022 2022 2022 ...
## $ definitivo2 : chr [1:157] "notifi" "notifi" "notifi" "notifi" ...
## $ fcaso : Factor w/ 2 levels "Control","Caso": 1 1 1 1 1 1 1 1 1 2 ...5 Centrality measures for numeric variables
v %>% group_by(ano) %>%
dplyr::select(Nuevo_dias_entre_notifi_visita,
Nuevo_dias_entre_enf_notifi,
poblacion,
Log_poblacion) %>%
descr()## Adding missing grouping variables: `ano`##
## ## Basic descriptive statistics
##
##
## Grouped by: 2022
##
## var type label n
## Nuevo_dias_entre_notifi_visita numeric Nuevo_dias_entre_notifi_visita 30
## Nuevo_dias_entre_enf_notifi numeric Nuevo_dias_entre_enf_notifi 30
## poblacion numeric poblacion 30
## NA.prc mean sd se md trimmed range iqr skew
## 0 0.43 0.50 0.09 0 0.42 1 (0-1) 1.00 0.28
## 0 17.43 23.10 4.22 8 12.25 101 (1-102) 16.25 2.40
## 0 15579.93 60570.68 11058.64 45 1972.88 329999 (1-330000) 751.50 5.16
##
##
## Grouped by: 2023
##
## var type label n
## Nuevo_dias_entre_notifi_visita numeric Nuevo_dias_entre_notifi_visita 109
## Nuevo_dias_entre_enf_notifi numeric Nuevo_dias_entre_enf_notifi 109
## poblacion numeric poblacion 109
## NA.prc mean sd se md trimmed range iqr skew
## 0 0.59 1.55 0.15 0 0.26 12 (0-12) 1 5.14
## 0 5.96 8.03 0.77 4 4.56 62 (0-62) 5 3.95
## 0 11848.33 57383.74 5496.37 20 1095.07 536910 (0-536910) 145 7.71
##
##
## Grouped by: 2024
##
## var type label n
## Nuevo_dias_entre_notifi_visita numeric Nuevo_dias_entre_notifi_visita 18
## Nuevo_dias_entre_enf_notifi numeric Nuevo_dias_entre_enf_notifi 18
## poblacion numeric poblacion 18
## NA.prc mean sd se md trimmed range iqr skew
## 0 0.28 0.46 0.11 0.0 0.25 1 (0-1) 0.75 1.08
## 0 10.78 17.39 4.10 4.5 8.00 76 (-5-71) 9.25 2.84
## 0 13801.67 48919.50 11530.44 50.5 2568.75 207330 (0-207330) 463.00 4.086 Tipo de explotacion
table(v$Nuevo_tipo, v$caso_o_notifi)##
## 0 1
## ENGORDE 17 0
## PONEDORA 12 15
## SILVESTRES 2 0
## TRASPATIO 93 18v$Nuevo_tipo <- gsub("SILVESTRES", "TRASPATIO", v$Nuevo_tipo)
v$ftipo <- as.factor(v$Nuevo_tipo)
tabela <- xtabs(~ v$ftipo + v$fcaso)
addmargins(tabela)## v$fcaso
## v$ftipo Control Caso Sum
## ENGORDE 17 0 17
## PONEDORA 12 15 27
## TRASPATIO 95 18 113
## Sum 124 33 157round(prop.table(tabela, 1),2)## v$fcaso
## v$ftipo Control Caso
## ENGORDE 1.00 0.00
## PONEDORA 0.44 0.56
## TRASPATIO 0.84 0.16fit.tipo <- glm(v$fcaso ~ v$ftipo, family = binomial(logit))
summary(fit.tipo)##
## Call:
## glm(formula = v$fcaso ~ v$ftipo, family = binomial(logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -18.57 1581.97 -0.012 0.991
## v$ftipoPONEDORA 18.79 1581.97 0.012 0.991
## v$ftipoTRASPATIO 16.90 1581.97 0.011 0.991
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 161.46 on 156 degrees of freedom
## Residual deviance: 136.20 on 154 degrees of freedom
## AIC: 142.2
##
## Number of Fisher Scoring iterations: 17cc(v$fcaso, v$ftipo)##
## v$ftipo
## v$fcaso ENGORDE PONEDORA TRASPATIO
## Control 17 12 95
## Caso 0 15 18
##
## Odds ratio 1 Inf Inf
## lower 95% CI 3.9 0.7
## upper 95% CI Inf Inf## Warning in chisq.test(cctable): Chi-squared approximation may be incorrect
## Warning in chisq.test(cctable): Chi-squared approximation may be incorrect## Warning in chisq.test(cctable, correct = FALSE): Chi-squared approximation may
## be incorrect## Chi-squared = 25.687 , 2 d.f., P value = 0
## Fisher's exact test (2-sided) P value = 0
##
## Cell counts too small - graph not shown
## chisq.test(tabela, correct = F)## Warning in chisq.test(tabela, correct = F): Chi-squared approximation may be
## incorrect##
## Pearson's Chi-squared test
##
## data: tabela
## X-squared = 25.687, df = 2, p-value = 2.644e-06ftipo <- v %>% group_by(fcaso, ftipo) %>%
summarise(Number_of_events=n()) %>%
ggplot() +
geom_col(aes(ftipo, Number_of_events, fill=fcaso)) +
scale_fill_brewer(palette = "Set2") +
# scale_fill_viridis_d(option = "C", direction = -1) +
xlab("Tipo de produccion") +
ylab("N eventos sanitarios")## `summarise()` has grouped output by 'fcaso'. You can override using the
## `.groups` argument.cc(v$fcaso, v$Nuevo_tipo)##
## v$Nuevo_tipo
## v$fcaso ENGORDE PONEDORA TRASPATIO
## Control 17 12 95
## Caso 0 15 18
##
## Odds ratio 1 Inf Inf
## lower 95% CI 3.9 0.7
## upper 95% CI Inf Inf## Warning in chisq.test(cctable): Chi-squared approximation may be incorrect## Warning in chisq.test(cctable): Chi-squared approximation may be incorrect## Warning in chisq.test(cctable, correct = FALSE): Chi-squared approximation may
## be incorrect## Chi-squared = 25.687 , 2 d.f., P value = 0
## Fisher's exact test (2-sided) P value = 0
##
## Cell counts too small - graph not shown
## 7 Bioseguridad
table(v$Nuevo_tipo, v$Nuevo_bioseguridad)##
## no si
## ENGORDE 11 6
## PONEDORA 13 14
## TRASPATIO 106 7v$fbio <- as.factor(v$Nuevo_bioseguridad)
tabela <- xtabs(~ v$fbio + v$fcaso)
addmargins(tabela)## v$fcaso
## v$fbio Control Caso Sum
## no 104 26 130
## si 20 7 27
## Sum 124 33 157round(prop.table(tabela, 1),2)## v$fcaso
## v$fbio Control Caso
## no 0.80 0.20
## si 0.74 0.26fit.bio <- glm(v$fcaso ~ v$fbio, family = binomial(logit))
summary(fit.bio)##
## Call:
## glm(formula = v$fcaso ~ v$fbio, family = binomial(logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3863 0.2193 -6.322 2.57e-10 ***
## v$fbiosi 0.3365 0.4908 0.685 0.493
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 161.46 on 156 degrees of freedom
## Residual deviance: 161.01 on 155 degrees of freedom
## AIC: 165.01
##
## Number of Fisher Scoring iterations: 4cc(v$fcaso, v$fbio)
##
## v$fbio
## v$fcaso no si Total
## Control 104 20 124
## Caso 26 7 33
## Total 130 27 157
##
## OR = 1.4
## 95% CI = 0.53, 3.66
## Chi-squared = 0.47, 1 d.f., P value = 0.492
## Fisher's exact test (2-sided) P value = 0.603exp(coef(fit.bio))## (Intercept) v$fbiosi
## 0.25 1.40chisq.test(tabela, correct = F)##
## Pearson's Chi-squared test
##
## data: tabela
## X-squared = 0.47292, df = 1, p-value = 0.4916fbio <- v %>% group_by(fcaso, fbio) %>%
summarise(Number_of_events=n()) %>%
ggplot() +
geom_col(aes(fbio, Number_of_events, fill=fcaso)) +
scale_fill_brewer(palette = "Set2") +
# scale_fill_viridis_d(option = "C", direction = -1) +
xlab("Bioseguridad") +
ylab("N eventos sanitarios")## `summarise()` has grouped output by 'fcaso'. You can override using the
## `.groups` argument.8 Notificador
table(v$fcaso, v$notificador)##
## AGROCALIDAD OTROS PROPIETARIO
## Control 11 20 93
## Caso 7 7 19v$fnot <- as.factor(v$notificador)
# Reordering reference value
v$fnot <- gsub("PROPIETARIO", "1PROPIETARIO", v$fnot)
tabela <- xtabs(~ v$fnot + v$fcaso)
addmargins(tabela)## v$fcaso
## v$fnot Control Caso Sum
## 1PROPIETARIO 93 19 112
## AGROCALIDAD 11 7 18
## OTROS 20 7 27
## Sum 124 33 157round(prop.table(tabela, 1),2)## v$fcaso
## v$fnot Control Caso
## 1PROPIETARIO 0.83 0.17
## AGROCALIDAD 0.61 0.39
## OTROS 0.74 0.26fit.not <- glm(v$fcaso ~ v$fnot, family = binomial(logit))
summary(fit.not)##
## Call:
## glm(formula = v$fcaso ~ v$fnot, family = binomial(logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.5882 0.2518 -6.308 2.82e-10 ***
## v$fnotAGROCALIDAD 1.1362 0.5451 2.084 0.0371 *
## v$fnotOTROS 0.5383 0.5062 1.063 0.2876
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 161.46 on 156 degrees of freedom
## Residual deviance: 156.95 on 154 degrees of freedom
## AIC: 162.95
##
## Number of Fisher Scoring iterations: 4cc(v$fcaso, v$fnot)##
## v$fnot
## v$fcaso 1PROPIETARIO AGROCALIDAD OTROS
## Control 93 11 20
## Caso 19 7 7
##
## Odds ratio 1 3.08 1.71
## lower 95% CI 0.89 0.53
## upper 95% CI 10.1 4.99## Warning in chisq.test(cctable): Chi-squared approximation may be incorrect
## Warning in chisq.test(cctable): Chi-squared approximation may be incorrect## Warning in chisq.test(cctable, correct = FALSE): Chi-squared approximation may
## be incorrect## Chi-squared = 4.963 , 2 d.f., P value = 0.084
## Fisher's exact test (2-sided) P value = 0.085
exp(coef(fit.not))## (Intercept) v$fnotAGROCALIDAD v$fnotOTROS
## 0.2043011 3.1148325 1.7131579chisq.test(tabela, correct = F)## Warning in chisq.test(tabela, correct = F): Chi-squared approximation may be
## incorrect##
## Pearson's Chi-squared test
##
## data: tabela
## X-squared = 4.9632, df = 2, p-value = 0.08361fnot <- v %>% group_by(fcaso, fnot) %>%
summarise(Number_of_events=n()) %>%
ggplot() +
geom_col(aes(fnot, Number_of_events, fill=fcaso)) +
scale_fill_brewer(palette = "Set2") +
# scale_fill_viridis_d(option = "C", direction = -1) +
xlab("Notificador") +
ylab("N eventos sanitarios")## `summarise()` has grouped output by 'fcaso'. You can override using the
## `.groups` argument.9 Sindrome presuntivo
table(v$fcaso, v$síndrome_presuntivo)##
## De Cabeza Inchada De Muerte Súbita Gastroentérico Respiratorio Agudo
## Control 1 13 8 72
## Caso 0 9 0 18
##
## Respiratorio Crónico
## Control 30
## Caso 6v$fsin <- as.factor(v$síndrome_presuntivo)
# Reordering reference value
v$fsin <- gsub("Respiratorio Crónico", "Respiratorio", v$fsin)
v$fsin <- gsub("Respiratorio Agudo", "Respiratorio", v$fsin)
v$fsin <- gsub("De Cabeza Inchada", "De Muerte Súbita", v$fsin)
v$fsin <- gsub("Respiratorio", "1Respiratorio", v$fsin)
tabela <- xtabs(~ v$fsin + v$fcaso)
addmargins(tabela)## v$fcaso
## v$fsin Control Caso Sum
## 1Respiratorio 102 24 126
## De Muerte Súbita 14 9 23
## Gastroentérico 8 0 8
## Sum 124 33 157round(prop.table(tabela, 1),2)## v$fcaso
## v$fsin Control Caso
## 1Respiratorio 0.81 0.19
## De Muerte Súbita 0.61 0.39
## Gastroentérico 1.00 0.00fit.sin <- glm(v$fcaso ~ v$fsin, family = binomial(logit))
summary(fit.sin)##
## Call:
## glm(formula = v$fcaso ~ v$fsin, family = binomial(logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.4469 0.2269 -6.378 1.8e-10 ***
## v$fsinDe Muerte Súbita 1.0051 0.4837 2.078 0.0377 *
## v$fsinGastroentérico -16.1191 1398.7210 -0.012 0.9908
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 161.46 on 156 degrees of freedom
## Residual deviance: 153.49 on 154 degrees of freedom
## AIC: 159.49
##
## Number of Fisher Scoring iterations: 16cc(v$fcaso, v$fsin)##
## v$fsin
## v$fcaso 1Respiratorio De Muerte Súbita Gastroentérico
## Control 102 14 8
## Caso 24 9 0
##
## Odds ratio 1 2.71 0
## lower 95% CI 0.92 0
## upper 95% CI 7.69 2.69## Warning in chisq.test(cctable): Chi-squared approximation may be incorrect
## Warning in chisq.test(cctable): Chi-squared approximation may be incorrect## Warning in chisq.test(cctable, correct = FALSE): Chi-squared approximation may
## be incorrect## Chi-squared = 6.969 , 2 d.f., P value = 0.031
## Fisher's exact test (2-sided) P value = 0.032
##
## Cell counts too small - graph not shown
## exp(coef(fit.sin))## (Intercept) v$fsinDe Muerte Súbita v$fsinGastroentérico
## 2.352941e-01 2.732143e+00 9.989467e-08chisq.test(tabela, correct = F)## Warning in chisq.test(tabela, correct = F): Chi-squared approximation may be
## incorrect##
## Pearson's Chi-squared test
##
## data: tabela
## X-squared = 6.9686, df = 2, p-value = 0.03068fsin <- v %>% group_by(fcaso, fsin) %>%
summarise(Number_of_events=n()) %>%
ggplot() +
geom_col(aes(fsin, Number_of_events, fill=fcaso)) +
scale_fill_brewer(palette = "Set2") +
# scale_fill_viridis_d(option = "C", direction = -1) +
xlab("Sindrome presuntivo") +
ylab("N eventos sanitarios")## `summarise()` has grouped output by 'fcaso'. You can override using the
## `.groups` argument.10 Tipo production
table(v$fcaso, v$Nuevo_tipo_comercio)##
## COMERCIAL >500 SILVESTRE TRASPATIO <500
## Control 22 2 100
## Caso 15 0 18v$fprod <- as.factor(v$Nuevo_tipo_comercio)
tabela <- xtabs(~ v$fprod + v$fcaso)
addmargins(tabela)## v$fcaso
## v$fprod Control Caso Sum
## COMERCIAL >500 22 15 37
## SILVESTRE 2 0 2
## TRASPATIO <500 100 18 118
## Sum 124 33 157round(prop.table(tabela, 1),2)## v$fcaso
## v$fprod Control Caso
## COMERCIAL >500 0.59 0.41
## SILVESTRE 1.00 0.00
## TRASPATIO <500 0.85 0.15fit.not <- glm(v$fcaso ~ v$fprod, family = binomial(logit))
summary(fit.not)##
## Call:
## glm(formula = v$fcaso ~ v$fprod, family = binomial(logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3830 0.3348 -1.144 0.25271
## v$fprodSILVESTRE -15.1831 1029.1215 -0.015 0.98823
## v$fprodTRASPATIO <500 -1.3318 0.4215 -3.160 0.00158 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 161.46 on 156 degrees of freedom
## Residual deviance: 150.75 on 154 degrees of freedom
## AIC: 156.75
##
## Number of Fisher Scoring iterations: 14cc(v$fcaso, v$fprod)##
## v$fprod
## v$fcaso COMERCIAL >500 SILVESTRE TRASPATIO <500
## Control 22 2 100
## Caso 15 0 18
##
## Odds ratio 1 0 0.27
## lower 95% CI 0 0.11
## upper 95% CI 8.55 0.66## Warning in chisq.test(cctable): Chi-squared approximation may be incorrect
## Warning in chisq.test(cctable): Chi-squared approximation may be incorrect## Warning in chisq.test(cctable, correct = FALSE): Chi-squared approximation may
## be incorrect## Chi-squared = 11.388 , 2 d.f., P value = 0.003
## Fisher's exact test (2-sided) P value = 0.003
##
## Cell counts too small - graph not shown
## exp(coef(fit.not))## (Intercept) v$fprodSILVESTRE v$fprodTRASPATIO <500
## 6.818182e-01 2.547264e-07 2.640000e-01chisq.test(tabela, correct = F)## Warning in chisq.test(tabela, correct = F): Chi-squared approximation may be
## incorrect##
## Pearson's Chi-squared test
##
## data: tabela
## X-squared = 11.388, df = 2, p-value = 0.003366fprod <- v %>% group_by(fcaso, fprod) %>%
summarise(Number_of_events=n()) %>%
ggplot() +
geom_col(aes(fprod, Number_of_events, fill=fcaso)) +
scale_fill_brewer(palette = "Set2") +
# scale_fill_viridis_d(option = "C", direction = -1) +
xlab("Tipo de produccion") +
ylab("N eventos sanitarios")## `summarise()` has grouped output by 'fcaso'. You can override using the
## `.groups` argument.ymd(v$f_1_ave_enferma) - ymd(v$f_notificación)## Time differences in days
## [1] -51 -11 -4 -16 -7 -6 -5 -5 -102 -11 -13 -1 -30 -5 -6
## [16] -23 -12 -1 -3 -8 -12 -5 -4 -3 -8 -2 -4 -23 -23 -2
## [31] -3 -1 -15 -2 -1 0 -10 -71 -3 -2 -7 0 -4 -4 -2
## [46] -14 -2 -3 -20 -1 -5 -11 -14 -4 -1 -52 -6 -1 -4 -9
## [61] -5 -3 -2 -2 -3 -28 -4 -2 -8 -6 -3 -4 -32 -5 -9
## [76] -11 -15 -3 -4 0 -1 -1 -5 -2 -14 -7 -1 -16 0 -2
## [91] -27 -3 0 -6 -7 -4 -15 0 0 0 -7 -3 -14 -13 -9
## [106] -4 -2 0 0 0 -12 -11 -6 -5 -5 -5 -1 -5 -1 -62
## [121] 0 0 -7 -4 -6 -7 -2 -2 -2 -3 -3 -7 0 0 -8
## [136] -3 -1 -10 -5 -4 -3 -71 0 -1 -9 -5 -2 -3 -10 -14
## [151] -13 -4 5 -13 -3 -36 -811 Resume of variables
library(ggpubr)## Warning: package 'ggpubr' was built under R version 4.3.3arr <- ggarrange(fbio,
fnot,
fsin,
ftipo,
fprod,
nrow = 3,
ncol = 2,
common.legend = TRUE,
font.label = list(size = 6, color = "black", family = NULL),
font.text = list(size = 6, color = "black", family = NULL),
legend = "bottom")
arr
11.1 R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.0011.2 Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.