Here is some example text
# colheita --> entrega
# Time from sample collection to receiving PCR result
#
# -- Some facilities (e.g. Munhava) have PCR facilities on-site
# Structural zeros in the distribution of time-to-event --> convergence issues
# -- Check out zero-inflated negative binomial random effects regression
# Link: https://stats.stackexchange.com/questions/38195/zero-inflated-negative-binomial-mixed-effects-model-in-r
#
time_fit1 <- glmer(
dias_colheita_entrega_pcr1 ~ intervencao + coorte + idade + (1|US),
data = filter(df2, dias_colheita_entrega_pcr1 > 3), # do this with zero-inflated model
family = poisson(link = "log")
)
summary(time_fit1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: dias_colheita_entrega_pcr1 ~ intervencao + coorte + idade + (1 |
## US)
## Data: filter(df2, dias_colheita_entrega_pcr1 > 3)
##
## AIC BIC logLik deviance df.resid
## 4421.1 4439.2 -2203.6 4407.1 91
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -11.3828 -4.5269 -0.9849 3.5893 21.8875
##
## Random effects:
## Groups Name Variance Std.Dev.
## US (Intercept) 0.09371 0.3061
## Number of obs: 98, groups: US, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.078964 0.132639 30.752 < 2e-16 ***
## intervencao 0.029260 0.040096 0.730 0.466
## coorte1a 0.537288 0.053783 9.990 < 2e-16 ***
## coorte2a 0.307123 0.052455 5.855 4.77e-09 ***
## coorte3a -0.034512 0.060238 -0.573 0.567
## idade 0.055505 0.003002 18.490 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intrvn coort1 coort2 coort3
## intervencao 0.027
## coorte1a -0.247 -0.425
## coorte2a -0.233 -0.567 0.815
## coorte3a -0.226 -0.679 0.790 0.891
## idade -0.142 -0.068 0.083 -0.064 0.094# negative binomial version doesn't converge
#
# time_fit2 <- glmer.nb(
# # formula = dias_colheita_entrega_pcr1 ~ intervencao + coorte + (1|US),
# formula = dias_colheita_entrega_pcr1 ~ intervencao + coorte + idade + (1|US),
# data = filter(df2, dias_colheita_entrega_pcr1 > 3 & dias_colheita_entrega_pcr1 < 180)
# )
# summary(time_fit2)
# get_ci(time_fit2, variable = "intervencao", exponentiate = TRUE)
# nascimento --> cpp
# Time from birth to starting CPP
time_fit3 <- glmer(
dias_nascimento_cpp ~ intervencao + coorte + idade + (1|US),
# data = df2,
data = filter(df2, dias_nascimento_cpp >= 0 & dias_nascimento_cpp < 90),
family = poisson(link = "log")
)
summary(time_fit3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: dias_nascimento_cpp ~ intervencao + coorte + idade + (1 | US)
## Data: filter(df2, dias_nascimento_cpp >= 0 & dias_nascimento_cpp <
## 90)
##
## AIC BIC logLik deviance df.resid
## 448.0 461.7 -217.0 434.0 45
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7395 -1.3844 -0.6403 0.3901 8.6110
##
## Random effects:
## Groups Name Variance Std.Dev.
## US (Intercept) 0.3539 0.5949
## Number of obs: 52, groups: US, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.50908 0.33463 7.498 6.47e-14 ***
## intervencao 0.64192 0.30054 2.136 0.032687 *
## coorte1a -0.77297 0.36766 -2.102 0.035518 *
## coorte2a -0.84358 0.34054 -2.477 0.013243 *
## coorte3a -1.37845 0.39145 -3.521 0.000429 ***
## idade 0.02021 0.01471 1.374 0.169571
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intrvn coort1 coort2 coort3
## intervencao 0.188
## coorte1a -0.476 -0.730
## coorte2a -0.491 -0.752 0.897
## coorte3a -0.479 -0.817 0.896 0.920
## idade -0.293 0.042 -0.023 -0.003 0.016time_fit4 <- glmer.nb(
formula = dias_nascimento_cpp ~ intervencao + coorte + idade + (1|US),
data = filter(df2, dias_nascimento_cpp >= 0 & dias_nascimento_cpp < 90)
)
summary(time_fit4)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(1.5162) ( log )
## Formula: dias_nascimento_cpp ~ intervencao + coorte + idade + (1 | US)
## Data: filter(df2, dias_nascimento_cpp >= 0 & dias_nascimento_cpp <
## 90)
##
## AIC BIC logLik deviance df.resid
## 326.9 342.5 -155.4 310.9 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1159 -0.7177 -0.3057 0.0973 3.8341
##
## Random effects:
## Groups Name Variance Std.Dev.
## US (Intercept) 0.1859 0.4311
## Number of obs: 52, groups: US, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.12365 0.49948 4.252 2.12e-05 ***
## intervencao 0.26771 0.54206 0.494 0.621
## coorte1a -0.30773 0.71415 -0.431 0.667
## coorte2a -0.31986 0.63446 -0.504 0.614
## coorte3a -0.81226 0.73678 -1.102 0.270
## idade 0.03190 0.03672 0.869 0.385
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intrvn coort1 coort2 coort3
## intervencao 0.164
## coorte1a -0.564 -0.576
## coorte2a -0.613 -0.660 0.775
## coorte3a -0.594 -0.788 0.796 0.863
## idade -0.403 -0.020 -0.088 -0.012 0.022get_ci(time_fit4, variable = "intervencao", exponentiate = TRUE)##
## Exponentiated 'intervencao': 1.31 (95% CI: 0.452, 3.78)# cpp --> ccr
# Time from starting CPP to starting CCR
time_fit5 <- glmer(
dias_cpp_ccr ~ intervencao + coorte + idade + (1|US),
data = filter(df2, dias_cpp_ccr >= 0 & dias_cpp_ccr < 60),
family = poisson(link = "log")
)
summary(time_fit5)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: dias_cpp_ccr ~ intervencao + coorte + idade + (1 | US)
## Data: filter(df2, dias_cpp_ccr >= 0 & dias_cpp_ccr < 60)
##
## AIC BIC logLik deviance df.resid
## 341.7 353.2 -163.8 327.7 31
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9046 -0.8258 0.2614 1.2140 2.3882
##
## Random effects:
## Groups Name Variance Std.Dev.
## US (Intercept) 0.06138 0.2478
## Number of obs: 38, groups: US, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.498424 0.181787 19.245 < 2e-16 ***
## intervencao -0.572831 0.147176 -3.892 9.94e-05 ***
## coorte1a -0.067896 0.179543 -0.378 0.7053
## coorte2a 0.001714 0.175566 0.010 0.9922
## coorte3a 0.379089 0.201266 1.884 0.0596 .
## idade -0.023735 0.010124 -2.344 0.0191 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intrvn coort1 coort2 coort3
## intervencao 0.027
## coorte1a -0.472 -0.377
## coorte2a -0.526 -0.542 0.636
## coorte3a -0.539 -0.701 0.683 0.848
## idade -0.372 -0.030 0.123 0.008 0.202time_fit6 <- glmer.nb(
formula = dias_cpp_ccr ~ intervencao + coorte + idade + (1|US),
data = filter(df2, dias_cpp_ccr >= 0 & dias_cpp_ccr < 60)
)
summary(time_fit6)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(6.2495) ( log )
## Formula: dias_cpp_ccr ~ intervencao + coorte + idade + (1 | US)
## Data: filter(df2, dias_cpp_ccr >= 0 & dias_cpp_ccr < 60)
##
## AIC BIC logLik deviance df.resid
## 315.2 328.3 -149.6 299.2 30
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0481 -0.3475 0.1855 0.6161 1.4969
##
## Random effects:
## Groups Name Variance Std.Dev.
## US (Intercept) 0.05545 0.2355
## Number of obs: 38, groups: US, 5
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.59146 0.34064 10.543 <2e-16 ***
## intervencao -0.66012 0.28894 -2.285 0.0223 *
## coorte1a -0.13879 0.40825 -0.340 0.7339
## coorte2a 0.03585 0.37876 0.095 0.9246
## coorte3a 0.44010 0.43151 1.020 0.3078
## idade -0.03738 0.02192 -1.706 0.0880 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intrvn coort1 coort2 coort3
## intervencao 0.015
## coorte1a -0.594 -0.348
## coorte2a -0.638 -0.528 0.663
## coorte3a -0.671 -0.658 0.692 0.865
## idade -0.411 0.049 0.093 -0.081 0.127get_ci(time_fit6, variable = "intervencao", exponentiate = TRUE)##
## Exponentiated 'intervencao': 0.517 (95% CI: 0.293, 0.91)# outcome: getting PCR result and start TARV on the same day
# -- exclude TDR
# Result: Compared to the control group, children in CCR receiving the intervention were
# more likely to start TARV on the same day that they received their
# PCR result (OR = 4.65 [95% CI: 1.48, 14.6]).
tmp <- df2[,
c("PCR", "TDR", "idade", "idade_start_TARV", "Data_Nascimento", "Data_Inscriacao_CCR", "Data_TDR",
"Data_1a_colheita_PCR", "Data_entrega_1o_PCR", "day_get_PCR",
"Data_inicio_TARV", "day_start_TARV", "same_day")]
df2_no_tdr <- filter(df2, TDR != 1)
df2_all_pcr <- filter(df2,
PCR == 1,
dias_colheita_entrega_pcr1 != 0,
birth_before_CCR,
CCR_before_PCRtake )
tmp2 <- tmp %>%
filter(TDR == 1)
# with age as covariate
time_fit7 <- glmer(
same_day ~ intervencao + coorte + idade + (1|US),
data = df2_all_pcr,
# data = df2_no_tdr,
family = binomial(link = "logit")
)
summary(time_fit7)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: same_day ~ intervencao + coorte + idade + (1 | US)
## Data: df2_all_pcr
##
## AIC BIC logLik deviance df.resid
## 101.7 117.4 -43.8 87.7 63
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5290 -1.0077 0.6705 0.7031 1.4618
##
## Random effects:
## Groups Name Variance Std.Dev.
## US (Intercept) 0 0
## Number of obs: 70, groups: US, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.042134 1.080558 0.039 0.969
## intervencao 1.524287 0.853167 1.787 0.074 .
## coorte1a -1.912857 1.473760 -1.298 0.194
## coorte2a -0.701937 1.239967 -0.566 0.571
## coorte3a -0.838917 1.359351 -0.617 0.537
## idade -0.007661 0.074423 -0.103 0.918
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intrvn coort1 coort2 coort3
## intervencao -0.032
## coorte1a -0.610 -0.430
## coorte2a -0.689 -0.453 0.743
## coorte3a -0.682 -0.623 0.767 0.869
## idade -0.379 0.084 -0.047 -0.151 0.002get_ci(time_fit7, variable = "intervencao", exponentiate = TRUE)##
## Exponentiated 'intervencao': 4.59 (95% CI: 0.863, 24.4)# without age as covariate
time_fit8 <- glmer(
same_day ~ intervencao + coorte + (1|US),
data = df2_no_tdr,
family = binomial(link = "logit")
)
summary(time_fit8)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: same_day ~ intervencao + coorte + (1 | US)
## Data: df2_no_tdr
##
## AIC BIC logLik deviance df.resid
## 199.9 218.1 -93.9 187.9 148
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3703 -1.0476 0.5233 0.7609 1.1841
##
## Random effects:
## Groups Name Variance Std.Dev.
## US (Intercept) 9e-16 3e-08
## Number of obs: 154, groups: US, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.8473 0.4879 1.736 0.08248 .
## intervencao 1.6331 0.5826 2.803 0.00506 **
## coorte1a -0.7543 0.7387 -1.021 0.30717
## coorte2a -1.1853 0.6930 -1.710 0.08719 .
## coorte3a -1.9338 0.8058 -2.400 0.01640 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intrvn coort1 coort2
## intervencao 0.000
## coorte1a -0.661 -0.149
## coorte2a -0.704 -0.560 0.549
## coorte3a -0.606 -0.723 0.508 0.831get_ci(time_fit8, variable = "intervencao", exponentiate = TRUE)##
## Exponentiated 'intervencao': 5.12 (95% CI: 1.63, 16)# stepped wedge analysis
# logistic regression, whether or not did fourth pick-up
fit1 <- glmer(
formula = lev2_dichot~intervencao + coorte + idade + (1|US),
# formula = lev4_dichot~intervencao + coorte + (1|US),
# formula = lev3_dichot~intervencao + idade + coorte + (1|US),
# formula = lev3_dichot~intervencao + coorte + (1|US),
# data = df,
data = df2,
family = binomial(link = "logit")
)
summary(fit1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: lev2_dichot ~ intervencao + coorte + idade + (1 | US)
## Data: df2
##
## AIC BIC logLik deviance df.resid
## 105.5 127.4 -45.8 91.5 161
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.5204 0.1863 0.2309 0.3355 0.9595
##
## Random effects:
## Groups Name Variance Std.Dev.
## US (Intercept) 0 0
## Number of obs: 168, groups: US, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.29997 0.72207 3.185 0.00145 **
## intervencao 0.78295 0.86212 0.908 0.36379
## coorte1a 0.43884 0.98276 0.447 0.65521
## coorte2a 0.80554 0.90587 0.889 0.37387
## coorte3a 0.41689 1.15954 0.360 0.71920
## idade -0.13859 0.06066 -2.285 0.02233 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intrvn coort1 coort2 coort3
## intervencao -0.040
## coorte1a -0.479 -0.184
## coorte2a -0.385 -0.536 0.466
## coorte3a -0.385 -0.741 0.428 0.708
## idade -0.603 0.067 0.019 -0.202 -0.018get_ci(fit1, variable = "intervencao", exponentiate = TRUE)##
## Exponentiated 'intervencao': 2.19 (95% CI: 0.404, 11.9)# Poisson regression (negative binomial version doesn't converge, funky distribution)
# -- Almost all have first few pick-ups, which is a good thing in general but not for getting a model to converge
fit2 <- glmer(
formula = lev_count~intervencao + idade + coorte + (1|US),
data = df,
family = poisson(link = "log"))
summary(fit2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: lev_count ~ intervencao + idade + coorte + (1 | US)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 586.5 608.3 -286.2 572.5 161
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.46020 0.06398 0.18734 0.26224 0.55086
##
## Random effects:
## Groups Name Variance Std.Dev.
## US (Intercept) 0 0
## Number of obs: 168, groups: US, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.288569 0.134311 9.594 <2e-16 ***
## intervencao 0.070956 0.132152 0.537 0.591
## idade -0.011052 0.009997 -1.106 0.269
## coorte1a 0.032154 0.176314 0.182 0.855
## coorte2a 0.005832 0.170677 0.034 0.973
## coorte3a -0.044762 0.190648 -0.235 0.814
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) intrvn idade coort1 coort2
## intervencao -0.043
## idade -0.430 0.099
## coorte1a -0.634 -0.240 0.031
## coorte2a -0.601 -0.586 -0.094 0.628
## coorte3a -0.572 -0.687 -0.005 0.606 0.852get_ci(fit2, variable = "intervencao", exponentiate = TRUE)##
## Exponentiated 'intervencao': 1.07 (95% CI: 0.829, 1.39)