CODES_growth_analysis1.R



CODES_growth_analysis1.Ruser2020-05-15rm(list = ls(all.names = TRUE))library(Amelia)## Loading required package: Rcpp## ## ## ## Amelia II: Multiple Imputation## ## (Version 1.7.6, built: 2019-11-24)## ## Copyright (C) 2005-2020 James Honaker, Gary King and Matthew Blackwell## ## Refer to for more information## ##library(bestNormalize)library(openxlsx)library(lubridate)## ## Attaching package: 'lubridate'## The following object is masked from 'package:base':## ## datelibrary(purrr)library(dplyr)## ## Attaching package: 'dplyr'## The following objects are masked from 'package:lubridate':## ## intersect, setdiff, union## The following objects are masked from 'package:stats':## ## filter, lag## The following objects are masked from 'package:base':## ## intersect, setdiff, setequal, unionlibrary(nlme)## ## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':## ## collapselibrary(rsq)library(lm.beta)library(tidyverse)## -- Attaching packages -------------------------------------------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.0 v readr 1.3.1## v tibble 2.1.3 v stringr 1.4.0## v tidyr 1.0.2 v forcats 0.4.0## -- Conflicts ----------------------------------------------------------------------------------- tidyverse_conflicts() --## x lubridate::as.difftime() masks base::as.difftime()## x nlme::collapse() masks dplyr::collapse()## x lubridate::date() masks base::date()## x dplyr::filter() masks stats::filter()## x lubridate::intersect() masks base::intersect()## x dplyr::lag() masks stats::lag()## x lubridate::setdiff() masks base::setdiff()## x lubridate::union() masks base::union()library(caret)## Loading required package: lattice## ## Attaching package: 'caret'## The following object is masked from 'package:purrr':## ## liftlibrary(leaps)library(ggplot2)library(reghelper)## ## Attaching package: 'reghelper'## The following object is masked from 'package:base':## ## betalibrary(ggpubr)## Loading required package: magrittr## ## Attaching package: 'magrittr'## The following object is masked from 'package:tidyr':## ## extract## The following object is masked from 'package:purrr':## ## set_nameslibrary(car)## Loading required package: carData## ## Attaching package: 'car'## The following object is masked from 'package:dplyr':## ## recode## The following object is masked from 'package:purrr':## ## somelibrary(missMDA)## Registered S3 methods overwritten by 'lme4':## method from## cooks.distance.influence.merMod car ## influence.merMod car ## dfbeta.influence.merMod car ## dfbetas.influence.merMod carlibrary(psych)## ## Attaching package: 'psych'## The following object is masked from 'package:car':## ## logit## The following object is masked from 'package:reghelper':## ## ICC## The following objects are masked from 'package:ggplot2':## ## %+%, alphasetwd("C:/Users/user/Dropbox/BCG-COVID-19/Ver 3/TF_SG")CV19 <- read.xlsx("Table_S1-9.xlsx", sheet=2)DATA <- CV19[, which (colnames(CV19) %in% c("Country_code","BCG.TB","m_cases","r_cases",#"popData2018.x",#"Positive.cases.per.test","temp","BCG.type","Current.BCG.vaccination","prevalence_overweight_2012","Fraction_urban_population","Population_ages_65_and_above_2018_percent","TB.incidents.per.10.thosand.mean_2000_2018","TB.Incidents.100"))]DATA_COMPLETE <- na.omit(DATA)dim(DATA_COMPLETE)## [1] 109 12#DATA_COMPLETE$pop2018_N <- predict(orderNorm(DATA_COMPLETE$popData2018.x))#hist(DATA_COMPLETE$pop2018_N)DATA_COMPLETE$TB_Z <- scale(DATA_COMPLETE$TB.incidents.per.10.thosand.mean_2000_2018)[,1]hist(DATA_COMPLETE$TB_Z)DATA_COMPLETE$Ob_N <- predict(orderNorm(DATA_COMPLETE$prevalence_overweight_2012))## Warning in orderNorm(DATA_COMPLETE$prevalence_overweight_2012): Ties in data, Normal distribution not guaranteedhist(DATA_COMPLETE$Ob_N)DATA_COMPLETE$temp_N <- predict(orderNorm(DATA_COMPLETE$temp))## Warning in orderNorm(DATA_COMPLETE$temp): Ties in data, Normal distribution not guaranteedhist(DATA_COMPLETE$temp_N)DATA_COMPLETE$frac_N <- predict(orderNorm(DATA_COMPLETE$Fraction_urban_population)) hist(DATA_COMPLETE$frac_N)DATA_COMPLETE$age_N <- predict(orderNorm(DATA_COMPLETE$Population_ages_65_and_above_2018_percent))hist(DATA_COMPLETE$age_N)DATA_COMPLETE$Yy <- ifelse(DATA_COMPLETE$BCG.TB == "Y_y", 1, 0)DATA_COMPLETE$Yn <- ifelse(DATA_COMPLETE$BCG.TB == "Y_n", 1, 0)DATA_COMPLETE$Nn <- ifelse(DATA_COMPLETE$BCG.TB == "N_n", 1, 0)###Analysis###shapiro.test(DATA_COMPLETE$m_cases)## ## Shapiro-Wilk normality test## ## data: DATA_COMPLETE$m_cases## W = 0.98452, p-value = 0.2394DATA_COMPLETE$m_cases_Z <- scale(DATA_COMPLETE$m_cases)[,1]summary(DATA_COMPLETE$r_cases)## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.5748 0.9171 0.9604 0.9351 0.9797 0.9966cases_gr_lm <- lm(m_cases ~ Current.BCG.vaccination + TB_Z +Ob_N + temp_N + age_N + frac_N,weights=DATA_COMPLETE$r_cases, data = DATA_COMPLETE)summary(cases_gr_lm, corr=T)## ## Call:## lm(formula = m_cases ~ Current.BCG.vaccination + TB_Z + Ob_N + ## temp_N + age_N + frac_N, data = DATA_COMPLETE, weights = DATA_COMPLETE$r_cases)## ## Weighted Residuals:## Min 1Q Median 3Q Max ## -0.081594 -0.017633 0.000579 0.015434 0.071780 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.1030093 0.0071712 14.364 < 2e-16 ***## Current.BCG.vaccinationY -0.0243129 0.0082116 -2.961 0.003817 ** ## TB_Z -0.0004480 0.0037284 -0.120 0.904593 ## Ob_N 0.0007051 0.0045013 0.157 0.875841 ## temp_N -0.0014532 0.0035765 -0.406 0.685351 ## age_N 0.0147254 0.0040473 3.638 0.000432 ***## frac_N 0.0011345 0.0044041 0.258 0.797240 ## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 0.02793 on 102 degrees of freedom## Multiple R-squared: 0.4031, Adjusted R-squared: 0.3679 ## F-statistic: 11.48 on 6 and 102 DF, p-value: 8.857e-10## ## Correlation of Coefficients:## (Intercept) Current.BCG.vaccinationY TB_Z Ob_N ## Current.BCG.vaccinationY -0.92 ## TB_Z 0.03 -0.03 ## Ob_N -0.04 0.05 0.40 ## temp_N 0.08 -0.09 -0.01 0.24## age_N -0.24 0.25 0.29 0.28## frac_N -0.13 0.14 -0.04 -0.61## temp_N age_N## Current.BCG.vaccinationY ## TB_Z ## Ob_N ## temp_N ## age_N 0.45 ## frac_N -0.13 -0.33shapiro.test(resid(cases_gr_lm)) #List of residuals## ## Shapiro-Wilk normality test## ## data: resid(cases_gr_lm)## W = 0.99563, p-value = 0.9825plot(density(resid(cases_gr_lm))) #A density plotrsq.partial(cases_gr_lm, adj=TRUE)## $adjustment## [1] TRUE## ## $variable## [1] "Current.BCG.vaccination" "TB_Z" ## [3] "Ob_N" "temp_N" ## [5] "age_N" "frac_N" ## ## $partial.rsq## [1] 0.070114013 -0.009661000 -0.009561084 -0.008172025 0.106196164## [6] -0.009147440lm.beta(cases_gr_lm)## ## Call:## lm(formula = m_cases ~ Current.BCG.vaccination + TB_Z + Ob_N + ## temp_N + age_N + frac_N, data = DATA_COMPLETE, weights = DATA_COMPLETE$r_cases)## ## Standardized Coefficients::## (Intercept) Current.BCG.vaccinationY TB_Z ## 0.00000000 -0.26069127 -0.01212401 ## Ob_N temp_N age_N ## 0.01905528 -0.03927589 0.39800398 ## frac_N ## 0.03066293Anova(cases_gr_lm)## Anova Table (Type II tests)## ## Response: m_cases## Sum Sq Df F value Pr(>F) ## Current.BCG.vaccination 0.006837 1 8.7663 0.0038167 ** ## TB_Z 0.000011 1 0.0144 0.9045925 ## Ob_N 0.000019 1 0.0245 0.8758412 ## temp_N 0.000129 1 0.1651 0.6853513 ## age_N 0.010325 1 13.2378 0.0004325 ***## frac_N 0.000052 1 0.0664 0.7972400 ## Residuals 0.079555 102 ## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1models_cases_gr <- regsubsets(m_cases ~ Current.BCG.vaccination + TB_Z +Ob_N + temp_N + age_N + frac_N,weights=DATA_COMPLETE$r_cases, data = DATA_COMPLETE, nvmax = 6)summary(models_cases_gr)## Subset selection object## Call: regsubsets.formula(m_cases ~ Current.BCG.vaccination + TB_Z + ## Ob_N + temp_N + age_N + frac_N, weights = DATA_COMPLETE$r_cases, ## data = DATA_COMPLETE, nvmax = 6)## 6 Variables (and intercept)## Forced in Forced out## Current.BCG.vaccinationY FALSE FALSE## TB_Z FALSE FALSE## Ob_N FALSE FALSE## temp_N FALSE FALSE## age_N FALSE FALSE## frac_N FALSE FALSE## 1 subsets of each size up to 6## Selection Algorithm: exhaustive## Current.BCG.vaccinationY TB_Z Ob_N temp_N age_N frac_N## 1 ( 1 ) " " " " " " " " "*" " " ## 2 ( 1 ) "*" " " " " " " "*" " " ## 3 ( 1 ) "*" " " "*" " " "*" " " ## 4 ( 1 ) "*" " " " " "*" "*" "*" ## 5 ( 1 ) "*" " " "*" "*" "*" "*" ## 6 ( 1 ) "*" "*" "*" "*" "*" "*"res.sum <- summary(models_cases_gr)data.frame( Adj.R2 = which.max(res.sum$adjr2), CP = which.min(res.sum$cp), BIC = which.min(res.sum$bic))## Adj.R2 CP BIC## 1 2 2 2coef(models_cases_gr, 1:5)## [[1]]## (Intercept) age_N ## 0.08340052 0.02123965 ## ## [[2]]## (Intercept) Current.BCG.vaccinationY age_N ## 0.10437798 -0.02602093 0.01616484 ## ## [[3]]## (Intercept) Current.BCG.vaccinationY Ob_N ## 0.103433425 -0.024847136 0.001888377 ## age_N ## 0.015821362 ## ## [[4]]## (Intercept) Current.BCG.vaccinationY temp_N ## 0.103134418 -0.024456369 -0.001672449 ## age_N frac_N ## 0.014704157 0.001754486 ## ## [[5]]## (Intercept) Current.BCG.vaccinationY Ob_N ## 0.1030371258 -0.0243413313 0.0009210781 ## temp_N age_N frac_N ## -0.0014565262 0.0148686862 0.0011135980vcov(models_cases_gr, 5)## (Intercept) Current.BCG.vaccinationY Ob_N## (Intercept) 5.088144e-05 -5.374896e-05 -1.781642e-06## Current.BCG.vaccinationY -5.374896e-05 6.673013e-05 2.106616e-06## Ob_N -1.781642e-06 2.106616e-06 1.686712e-05## temp_N 1.966519e-06 -2.497721e-06 3.954055e-06## age_N -7.218883e-06 8.563101e-06 3.012927e-06## frac_N -4.004879e-06 5.139757e-06 -1.173618e-05## temp_N age_N frac_N## (Intercept) 1.966519e-06 -7.218883e-06 -4.004879e-06## Current.BCG.vaccinationY -2.497721e-06 8.563101e-06 5.139757e-06## Ob_N 3.954055e-06 3.012927e-06 -1.173618e-05## temp_N 1.266826e-05 6.472628e-06 -2.106404e-06## age_N 6.472628e-06 1.481588e-05 -5.581528e-06## frac_N -2.106404e-06 -5.581528e-06 1.918090e-05cases_gr_lm_2 <- lm(m_cases ~ Current.BCG.vaccination + age_N,weights=DATA_COMPLETE$r_cases, data = DATA_COMPLETE)summary(cases_gr_lm_2, corr=T)## ## Call:## lm(formula = m_cases ~ Current.BCG.vaccination + age_N, data = DATA_COMPLETE, ## weights = DATA_COMPLETE$r_cases)## ## Weighted Residuals:## Min 1Q Median 3Q Max ## -0.082933 -0.019419 -0.000057 0.015442 0.072399 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.104378 0.006819 15.307 < 2e-16 ***## Current.BCG.vaccinationY -0.026021 0.007755 -3.355 0.0011 ** ## age_N 0.016165 0.003146 5.138 1.28e-06 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 0.02748 on 106 degrees of freedom## Multiple R-squared: 0.3995, Adjusted R-squared: 0.3882 ## F-statistic: 35.26 on 2 and 106 DF, p-value: 1.825e-12## ## Correlation of Coefficients:## (Intercept) Current.BCG.vaccinationY## Current.BCG.vaccinationY -0.92 ## age_N -0.45 0.48rsq.partial(cases_gr_lm_2, adj=TRUE)## $adjustment## [1] TRUE## ## $variable## [1] "Current.BCG.vaccination" "age_N" ## ## $partial.rsq## [1] 0.08749114 0.19180907lm.beta(cases_gr_lm_2)## ## Call:## lm(formula = m_cases ~ Current.BCG.vaccination + age_N, data = DATA_COMPLETE, ## weights = DATA_COMPLETE$r_cases)## ## Standardized Coefficients::## (Intercept) Current.BCG.vaccinationY age_N ## 0.0000000 -0.2790055 0.4369088DATA_d <- CV19[, which (colnames(CV19) %in% c("Country_code","BCG.TB","m_deaths","r_deaths",#"popData2018.x",#"Positive.cases.per.test","temp","BCG.type","Current.BCG.vaccination","prevalence_overweight_2012","Fraction_urban_population","Population_ages_65_and_above_2018_percent","TB.incidents.per.10.thosand.mean_2000_2018","TB.Incidents.100"))]DATA_COMPLETE_d <- na.omit(DATA_d)dim(DATA_COMPLETE_d)## [1] 66 12#DATA_COMPLETE_d$pop2018_N <- predict(orderNorm(DATA_COMPLETE_d$popData2018.x))#hist(DATA_COMPLETE_d$pop2018_N)DATA_COMPLETE_d$TB_Z <- scale(DATA_COMPLETE_d$TB.incidents.per.10.thosand.mean_2000_2018)[,1]hist(DATA_COMPLETE_d$TB_Z)DATA_COMPLETE_d$Ob_N <- predict(orderNorm(DATA_COMPLETE_d$prevalence_overweight_2012))## Warning in orderNorm(DATA_COMPLETE_d$prevalence_overweight_2012): Ties in data, Normal distribution not guaranteedhist(DATA_COMPLETE_d$Ob_N)DATA_COMPLETE_d$temp_N <- predict(orderNorm(DATA_COMPLETE_d$temp))## Warning in orderNorm(DATA_COMPLETE_d$temp): Ties in data, Normal distribution not guaranteedhist(DATA_COMPLETE_d$temp_N)DATA_COMPLETE_d$frac_N <- predict(orderNorm(DATA_COMPLETE_d$Fraction_urban_population)) hist(DATA_COMPLETE_d$frac_N)DATA_COMPLETE_d$age_N <- predict(orderNorm(DATA_COMPLETE_d$Population_ages_65_and_above_2018_percent))hist(DATA_COMPLETE_d$age_N)###Analysis###shapiro.test(DATA_COMPLETE_d$m_deaths)## ## Shapiro-Wilk normality test## ## data: DATA_COMPLETE_d$m_deaths## W = 0.71942, p-value = 6.339e-10DATA_COMPLETE_d$m_deaths_N <- predict(orderNorm(DATA_COMPLETE_d$m_deaths))DATA_COMPLETE_d$m_deaths_Z <- scale(DATA_COMPLETE_d$m_deaths)[,1]summary(DATA_COMPLETE_d$r_deaths)## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.7027 0.9095 0.9544 0.9375 0.9801 0.9939deaths_gr_lm <- lm(m_deaths_N~ Current.BCG.vaccination + TB_Z +Ob_N + temp_N + age_N + frac_N,weights=DATA_COMPLETE_d$r_deaths, data = DATA_COMPLETE_d)summary(deaths_gr_lm, corr=T)## ## Call:## lm(formula = m_deaths_N ~ Current.BCG.vaccination + TB_Z + Ob_N + ## temp_N + age_N + frac_N, data = DATA_COMPLETE_d, weights = DATA_COMPLETE_d$r_deaths)## ## Weighted Residuals:## Min 1Q Median 3Q Max ## -1.67101 -0.43334 -0.01121 0.51739 1.64505 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.501040 0.203256 2.465 0.01663 * ## Current.BCG.vaccinationY -0.715306 0.255391 -2.801 0.00688 **## TB_Z 0.037058 0.117829 0.315 0.75425 ## Ob_N 0.073942 0.125818 0.588 0.55898 ## temp_N -0.079472 0.128483 -0.619 0.53860 ## age_N 0.367053 0.134440 2.730 0.00833 **## frac_N 0.006404 0.126152 0.051 0.95969 ## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 0.7525 on 59 degrees of freedom## Multiple R-squared: 0.4316, Adjusted R-squared: 0.3738 ## F-statistic: 7.467 on 6 and 59 DF, p-value: 5.644e-06## ## Correlation of Coefficients:## (Intercept) Current.BCG.vaccinationY TB_Z Ob_N ## Current.BCG.vaccinationY -0.88 ## TB_Z 0.06 -0.07 ## Ob_N -0.01 0.01 0.33 ## temp_N 0.12 -0.14 -0.13 0.21## age_N -0.22 0.25 0.15 0.30## frac_N -0.22 0.26 0.12 -0.46## temp_N age_N## Current.BCG.vaccinationY ## TB_Z ## Ob_N ## temp_N ## age_N 0.52 ## frac_N -0.27 -0.22shapiro.test(resid(deaths_gr_lm)) #List of residuals## ## Shapiro-Wilk normality test## ## data: resid(deaths_gr_lm)## W = 0.99244, p-value = 0.9609plot(density(resid(deaths_gr_lm))) #A density plotrsq.partial(deaths_gr_lm, adj=TRUE)## $adjustment## [1] TRUE## ## $variable## [1] "Current.BCG.vaccination" "TB_Z" ## [3] "Ob_N" "temp_N" ## [5] "age_N" "frac_N" ## ## $partial.rsq## [1] 0.10239604 -0.01524709 -0.01103069 -0.01039705 0.09712290 -0.01690474lm.beta(deaths_gr_lm)## ## Call:## lm(formula = m_deaths_N ~ Current.BCG.vaccination + TB_Z + Ob_N + ## temp_N + age_N + frac_N, data = DATA_COMPLETE_d, weights = DATA_COMPLETE_d$r_deaths)## ## Standardized Coefficients::## (Intercept) Current.BCG.vaccinationY TB_Z ## 0.000000000 -0.331917234 0.037132336 ## Ob_N temp_N age_N ## 0.073928858 -0.079471559 0.367053357 ## frac_N ## 0.006403789Anova(deaths_gr_lm)## Anova Table (Type II tests)## ## Response: m_deaths_N## Sum Sq Df F value Pr(>F) ## Current.BCG.vaccination 4.442 1 7.8446 0.006881 **## TB_Z 0.056 1 0.0989 0.754247 ## Ob_N 0.196 1 0.3454 0.558983 ## temp_N 0.217 1 0.3826 0.538598 ## age_N 4.221 1 7.4542 0.008331 **## frac_N 0.001 1 0.0026 0.959686 ## Residuals 33.407 59 ## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1models_deaths_gr <- regsubsets(m_deaths_N ~ Current.BCG.vaccination + TB_Z +Ob_N + temp_N + age_N + frac_N,weights=DATA_COMPLETE_d$r_cases, data = DATA_COMPLETE_d, nvmax = 6)summary(models_deaths_gr)## Subset selection object## Call: regsubsets.formula(m_deaths_N ~ Current.BCG.vaccination + TB_Z + ## Ob_N + temp_N + age_N + frac_N, weights = DATA_COMPLETE_d$r_cases, ## data = DATA_COMPLETE_d, nvmax = 6)## 6 Variables (and intercept)## Forced in Forced out## Current.BCG.vaccinationY FALSE FALSE## TB_Z FALSE FALSE## Ob_N FALSE FALSE## temp_N FALSE FALSE## age_N FALSE FALSE## frac_N FALSE FALSE## 1 subsets of each size up to 6## Selection Algorithm: exhaustive## Current.BCG.vaccinationY TB_Z Ob_N temp_N age_N frac_N## 1 ( 1 ) " " " " " " " " "*" " " ## 2 ( 1 ) "*" " " " " " " "*" " " ## 3 ( 1 ) "*" " " "*" " " "*" " " ## 4 ( 1 ) "*" " " "*" "*" "*" " " ## 5 ( 1 ) "*" "*" "*" "*" "*" " " ## 6 ( 1 ) "*" "*" "*" "*" "*" "*"res.sum <- summary(models_deaths_gr)data.frame( Adj.R2 = which.max(res.sum$adjr2), CP = which.min(res.sum$cp), BIC = which.min(res.sum$bic))## Adj.R2 CP BIC## 1 2 2 2coef(models_deaths_gr, 1:5)## [[1]]## (Intercept) age_N ## 1.734723e-17 5.621013e-01 ## ## [[2]]## (Intercept) Current.BCG.vaccinationY age_N ## 0.5662117 -0.8123907 0.3867818 ## ## [[3]]## (Intercept) Current.BCG.vaccinationY Ob_N ## 0.53497164 -0.76757133 0.07571026 ## age_N ## 0.39432555 ## ## [[4]]## (Intercept) Current.BCG.vaccinationY Ob_N ## 0.52663508 -0.75560935 0.06488320 ## temp_N age_N ## -0.07569825 0.35214509 ## ## [[5]]## (Intercept) Current.BCG.vaccinationY TB_Z ## 0.53344617 -0.76538255 0.04321621 ## Ob_N temp_N age_N ## 0.08267392 -0.08021763 0.36113505vcov(models_deaths_gr, 5)## (Intercept) Current.BCG.vaccinationY TB_Z## (Intercept) 0.038411466 -0.041917950 0.002187231## Current.BCG.vaccinationY -0.041917950 0.060142968 -0.003138446## TB_Z 0.002187231 -0.003138446 0.013877920## Ob_N -0.002976471 0.004270026 0.005713093## temp_N 0.001446267 -0.002075219 -0.001451297## age_N -0.007378247 0.010585993 0.002886927## Ob_N temp_N age_N## (Intercept) -0.002976471 0.001446267 -0.007378247## Current.BCG.vaccinationY 0.004270026 -0.002075219 0.010585993## TB_Z 0.005713093 -0.001451297 0.002886927## Ob_N 0.012639267 0.001577942 0.003394652## temp_N 0.001577942 0.015361218 0.008173081## age_N 0.003394652 0.008173081 0.017384373deaths_gr_lm_2 <- lm(m_deaths_N ~ Current.BCG.vaccination + age_N,weights=DATA_COMPLETE_d$r_cases, data = DATA_COMPLETE_d)summary(deaths_gr_lm_2, corr=T)## ## Call:## lm(formula = m_deaths_N ~ Current.BCG.vaccination + age_N, data = DATA_COMPLETE_d, ## weights = DATA_COMPLETE_d$r_cases)## ## Residuals:## Min 1Q Median 3Q Max ## -2.0240 -0.4080 0.0123 0.5525 1.7430 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.5662 0.1873 3.023 0.003618 ** ## Current.BCG.vaccinationY -0.8124 0.2321 -3.500 0.000860 ***## age_N 0.3868 0.1077 3.591 0.000645 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 0.7671 on 63 degrees of freedom## Multiple R-squared: 0.4273, Adjusted R-squared: 0.4091 ## F-statistic: 23.51 on 2 and 63 DF, p-value: 2.367e-08## ## Correlation of Coefficients:## (Intercept) Current.BCG.vaccinationY## Current.BCG.vaccinationY -0.86 ## age_N -0.40 0.47rsq.partial(deaths_gr_lm_2, adj=TRUE)## $adjustment## [1] TRUE## ## $variable## [1] "Current.BCG.vaccination" "age_N" ## ## $partial.rsq## [1] 0.1495184 0.1567605lm.beta(deaths_gr_lm_2)## ## Call:## lm(formula = m_deaths_N ~ Current.BCG.vaccination + age_N, data = DATA_COMPLETE_d, ## weights = DATA_COMPLETE_d$r_cases)## ## Standardized Coefficients::## (Intercept) Current.BCG.vaccinationY age_N ## 0.0000000 -0.3769667 0.3867818######PCA#cases_PCA_lm <- lm(m_cases ~ Current.BCG.vaccination + TB_Z +RC1 + RC2 + RC3,weights=CV19$r_cases, data = CV19)summary(cases_PCA_lm, corr=T)## ## Call:## lm(formula = m_cases ~ Current.BCG.vaccination + TB_Z + RC1 + ## RC2 + RC3, data = CV19, weights = CV19$r_cases)## ## Weighted Residuals:## Min 1Q Median 3Q Max ## -0.059369 -0.017509 -0.000252 0.016160 0.066387 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.0975696 0.0073366 13.299 < 2e-16 ***## Current.BCG.vaccinationY -0.0186245 0.0082860 -2.248 0.0263 * ## TB_Z -0.0027056 0.0032983 -0.820 0.4136 ## RC1 0.0153722 0.0036100 4.258 4.01e-05 ***## RC2 0.0017034 0.0026717 0.638 0.5249 ## RC3 0.0008631 0.0026472 0.326 0.7449 ## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 0.02907 on 125 degrees of freedom## Multiple R-squared: 0.3523, Adjusted R-squared: 0.3264 ## F-statistic: 13.6 on 5 and 125 DF, p-value: 1.356e-10## ## Correlation of Coefficients:## (Intercept) Current.BCG.vaccinationY TB_Z RC1 RC2 ## Current.BCG.vaccinationY -0.93 ## TB_Z 0.03 -0.02 ## RC1 -0.45 0.48 0.47 ## RC2 -0.06 0.07 0.05 0.06 ## RC3 -0.12 0.12 -0.03 0.05 0.01shapiro.test(resid(cases_PCA_lm)) #List of residuals## ## Shapiro-Wilk normality test## ## data: resid(cases_PCA_lm)## W = 0.98397, p-value = 0.1267plot(density(resid(cases_PCA_lm))) #A density plotrsq.partial(cases_PCA_lm, adj=TRUE)## $adjustment## [1] TRUE## ## $variable## [1] "Current.BCG.vaccination" "TB_Z" ## [3] "RC1" "RC2" ## [5] "RC3" ## ## $partial.rsq## [1] 0.031158143 -0.002602732 0.119696173 -0.004732732 -0.007143555lm.beta(cases_PCA_lm)## ## Call:## lm(formula = m_cases ~ Current.BCG.vaccination + TB_Z + RC1 + ## RC2 + RC3, data = CV19, weights = CV19$r_cases)## ## Standardized Coefficients::## (Intercept) Current.BCG.vaccinationY TB_Z ## 0.00000000 -0.19059045 -0.07249564 ## RC1 RC2 RC3 ## 0.41189235 0.04564151 0.02312588kruskal.test(RC1 ~ Current.BCG.vaccination, data = CV19)## ## Kruskal-Wallis rank sum test## ## data: RC1 by Current.BCG.vaccination## Kruskal-Wallis chi-squared = 44.445, df = 1, p-value = 2.617e-11kruskal.test(RC2 ~ Current.BCG.vaccination, data = CV19)## ## Kruskal-Wallis rank sum test## ## data: RC2 by Current.BCG.vaccination## Kruskal-Wallis chi-squared = 0.70711, df = 1, p-value = 0.4004kruskal.test(RC3 ~ Current.BCG.vaccination, data = CV19)## ## Kruskal-Wallis rank sum test## ## data: RC3 by Current.BCG.vaccination## Kruskal-Wallis chi-squared = 0.76947, df = 1, p-value = 0.3804dataset_1 <- na.omit(CV19[, which (colnames(CV19) %in% c("m_deaths_N","r_deaths","Current.BCG.vaccination","TB_Z","RC1","RC2","RC3"))])dataset_1<-na.omit(dataset_1)deaths_PCA_lm <- lm(m_deaths_N ~ Current.BCG.vaccination + TB_Z +RC1 + RC2 + RC3,weights=dataset_1$r_deaths, data = dataset_1)summary(deaths_PCA_lm, corr=T)## ## Call:## lm(formula = m_deaths_N ~ Current.BCG.vaccination + TB_Z + RC1 + ## RC2 + RC3, data = dataset_1, weights = dataset_1$r_deaths)## ## Weighted Residuals:## Min 1Q Median 3Q Max ## -1.6699 -0.4195 -0.0299 0.4943 1.4742 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.03426 0.28204 0.121 0.90367 ## Current.BCG.vaccinationY -0.52897 0.26431 -2.001 0.04924 * ## TB_Z -0.16524 0.20407 -0.810 0.42086 ## RC1 0.50486 0.17467 2.890 0.00512 **## RC2 0.11484 0.10266 1.119 0.26710 ## RC3 0.12992 0.08920 1.457 0.14972 ## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Residual standard error: 0.7775 on 70 degrees of freedom## Multiple R-squared: 0.3681, Adjusted R-squared: 0.323 ## F-statistic: 8.155 on 5 and 70 DF, p-value: 4.148e-06## ## Correlation of Coefficients:## (Intercept) Current.BCG.vaccinationY TB_Z RC1 RC2 ## Current.BCG.vaccinationY -0.89 ## TB_Z 0.28 -0.11 ## RC1 -0.61 0.55 0.30 ## RC2 -0.13 0.09 -0.13 0.02 ## RC3 -0.27 0.21 -0.10 0.19 0.07shapiro.test(resid(deaths_PCA_lm)) #List of residuals## ## Shapiro-Wilk normality test## ## data: resid(deaths_PCA_lm)## W = 0.97316, p-value = 0.1068plot(density(resid(deaths_PCA_lm))) #A density plotrsq.partial(deaths_PCA_lm, adj=TRUE)## $adjustment## [1] TRUE## ## $variable## [1] "Current.BCG.vaccination" "TB_Z" ## [3] "RC1" "RC2" ## [5] "RC3" ## ## $partial.rsq## [1] 0.040607159 -0.004874348 0.093863152 0.003528800 0.015549178lm.beta(deaths_PCA_lm)## ## Call:## lm(formula = m_deaths_N ~ Current.BCG.vaccination + TB_Z + RC1 + ## RC2 + RC3, data = dataset_1, weights = dataset_1$r_deaths)## ## Standardized Coefficients::## (Intercept) Current.BCG.vaccinationY TB_Z ## 0.00000000 -0.24191645 -0.08513101 ## RC1 RC2 RC3 ## 0.36318781 0.10669203 0.14166647 ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download