Lưu trữ các biến giải thích

Tên của các biến giải thích sẽ được lưu trữ dưới dạng một véc tơ ký tự. Véc tơ này sẽ được đề cập về sau.

## define variables of interest 
explanatory_vars <- c("gender", "fever", "chills", "cough", "aches", "vomit")

Chuyển đổi sang số 1 và số 0

Sau đây, giá trị của các biến giải thích được chuyển đổi từ “có”/“không”, “nam”/“nữ” và “chết”/“sống” thành 1 / 0, để hợp với các đặc tính của mô hình hồi quy logistic. TĐể thực hiện việc này một cách hiệu quả, sử dụng hàm across() từ dplyr để chuyển đổi nhiều biến cùng một lúc. Để áp dụng cho mỗi biến, dùng hàm case_when() (cũng trong package dplyr) để chuyển đổi các giá trị cụ thể thành 1 và 0. Xem các mục về across() và case_when() trong chương Làm sạch số liệu và các hàm quan trọng).

Chú ý: dấu “.” bên dưới đại diện cho cột`````đang được xử lý trong hàmacross()` tại thời điểm đó.

## convert dichotomous variables to 0/1 
linelist <- linelist %>%  
  mutate(across(                                      
    .cols = all_of(c(explanatory_vars, "outcome")),  ## for each column listed and "outcome"
    .fns = ~case_when(                              
      . %in% c("m", "yes", "Death")   ~ 1,           ## recode male, yes and death to 1
      . %in% c("f", "no",  "Recover") ~ 0,           ## female, no and recover to 0
      TRUE                            ~ NA_real_)    ## otherwise set to missing
    )
  )

Loại bỏ các hàng có giá trị missing

Để bỏ các hàng có giá trị missing, dùng hàm drop_na() trong package tidyr. Tuy nhiên, chúng ta chỉ muốn thực hiện điều này cho các hàng có giá trị missing đối với các cột đang được quan tâm.

Trước hết, chúng ta phải đảm bảo rằng vectơ explanatory_vars bao gồm các biến age (age có thể tạo ra một lỗi trong thao tác của hàm case_when() trước đó, mà chỉ dành cho biến nhị phân). Sau đó chúng ta pipe bộ dữ liệu linelist tới hàm drop_na() để bỏ các hàng có giá trị missing cho biến outcome hoặc bất kỳ biển giải thích explanatory_vars nào.

Trước khi thực hiện các lệnh này, kiểm tra số hàng trong bộ số liệu linelist bằng hàm nrow(linelist).

## add in age_category to the explanatory vars 
explanatory_vars <- c(explanatory_vars, "age_cat")

## drop rows with missing information for variables of interest 
linelist <- linelist %>% 
  drop_na(any_of(c("outcome", explanatory_vars)))

Kiểm tra số hàng còn lại của linelist bằng hàm nrow(linelist).

Hồi quy tuyến tính

Hàm lm() trong base cho phép thực hiện hồi quy tuyến tính để đánh giá mối quan hệ giữa biến đầu ra dạng số (numeric) và các biến giải thích mà được giả định là có mối quan hệ tuyến tính.

Cung cấp phương trình dưới dạng công thức với tên của biến đầu ra và các biến giải thích được phân tách bằng dấu ngã ~. Bên cạnh đó, chỉ rõ bộ số liệu nào được sử dụng với data =. Kết quả của mô hình được định nghĩa dưới dạng đối tượng của R để sử dụng về sau.

lm_results <- lm(ht_cm ~ age, data = linelist)

Sau đó tóm tắt kết quả của mô hình bằng hàm summary() để xem các hệ số (ước tính), P-value, phần dư và các đo lường khác.

summary(lm_results)

## 
## Call:
## lm(formula = ht_cm ~ age, data = linelist)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -128.579  -15.854    1.177   15.887  175.483 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  69.9051     0.5979   116.9   <2e-16 ***
## age           3.4354     0.0293   117.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.75 on 4165 degrees of freedom
## Multiple R-squared:  0.7675, Adjusted R-squared:  0.7674 
## F-statistic: 1.375e+04 on 1 and 4165 DF,  p-value: < 2.2e-16

Ngoài ra, có thể dùng hàm tidy() trong package broom để xuất kết quả vào trong một bảng. Kết quả bên dưới cho chúng ta biết khi tăng thêm một tuổi thì chiều cao tăng 3,5 cm và mối quan hệ này có ý nghĩa thống kê.

tidy(lm_results)

## # A tibble: 2 x 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)    69.9     0.598       117.       0
## 2 age             3.44    0.0293      117.       0

Sau đó, có thể sử dụng kết quả hồi quy này để đưa vào ggplot. Để thực hiện điều này, trước tiên chúng ta đưa các giá trị quan sát và đường thẳng hồi quy (fitted line) vào một data frame bằng cách dùng hàm augment() trong package broom.

## pull the regression points and observed data in to one dataset
points <- augment(lm_results)

## plot the data using age as the x-axis 
ggplot(points, aes(x = age)) + 
  ## add points for height 
  geom_point(aes(y = ht_cm)) + 
  ## add your regression line 
  geom_line(aes(y = .fitted), colour = "red")

Bạn cũng có thể vẽ đường hồi quy tuyến tính đơn bằng package ggplot thông qua hàm geom_smooth().

## add your data to a plot 
 ggplot(linelist, aes(x = age, y = ht_cm)) + 
  ## show points
  geom_point() + 
  ## add a linear regression 
  geom_smooth(method = "lm", se = FALSE)

## `geom_smooth()` using formula 'y ~ x'

Xem thêm các hướng dẫn chi tiết trong mục Nguồn ở cuối chương này.

Hồi quy Logistic

Hàm glm() trong package stats (một phần của base R) được sử dụng để fit (chọn mô hình dự đoán tối ưu dựa trên số liệu quan sát) đối với Mô hình Tuyến tính Tổng quát (GLM).

glm() có thể được sử dụng cho cả hồi quy logistic đơn biến và đa biến (ví dụ như để tính tỷ số chênh OR). Sau đây là những thành phần chính của hàm:

# arguments for glm()
glm(formula, family, data, weights, subset, ...)

• formula = Mô hình được cung cấp cho glm() dưới dạng một phương trình với biến kết cục ở bên trái và biến giải thích ở bên phải dấu ngã ~.
  
• family = Xác định loại mô hình sẽ thực hiện. Đối với hồi quy logistic, sử dụng family = "binomial", đối với hồi quy poisson sử dụng family = "poisson". Các ví dụ khác được trình bày trong bảng bên dưới.
  
• data = Cụ thể bộ số liệu

Nếu cần, có thể cụ thể hàm liên kết bằng cú pháp family = familytype(link = "linkfunction")). Bạn có thể tìm đọc thêm về các họ hồi quy khác và các tùy chọn đối số như là weights = và subset = bằng cách gõ (?glm).

Khi thực hiện glm() , phổ biến nhất là lưu kết quả dưới dạng một đối tượng của R được đặt tên. Sau đó, có thể xuất kết quả ra console bằng cách sử dụng hàm summary() như được trình bày bên dưới, hoặc thực hiện các thao tác khác từ kết quả (ví dụ như lấy lũy thừa).

Nếu cần thực hiện một hồi quy nhị thức âm, có thể sử dụng package MASS. Hàn glm.nb() uses cũng sử dụng cùng cú pháp như glm(). Để xem qua các hồi quy khác, xem trên trang thống kê của UCLA.

Phân tích đơn biến sử dụng glm()

Trong ví dụ này, chúng tôi đánh giá mối liên quan giữa nhóm tuổi và biến kết cục tử vong (được mã hóa là 1 trong phần chuẩn bị). Bên dưới là một mô hình đơn biến của biến kết cục outcome theo age_cat. Chúng tôi lưu kết quả đầu ra được đặt tên là model và sau đó in kết quả đến console bằng hàm summary(). Lưu ý, các ước tính được tạo ra là các giá trị lôgarít của tỷ số chênh (log odds) và giá trị tham chiếu là giá trị đầu tiên của biến age_cat (“0-4”).

model <- glm(outcome ~ age_cat, family = "binomial", data = linelist)
summary(model)

## 
## Call:
## glm(formula = outcome ~ age_cat, family = "binomial", data = linelist)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.339  -1.278   1.024   1.080   1.354  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.233738   0.072805   3.210  0.00133 **
## age_cat5-9   -0.062898   0.101733  -0.618  0.53640   
## age_cat10-14  0.138204   0.107186   1.289  0.19726   
## age_cat15-19 -0.005565   0.113343  -0.049  0.96084   
## age_cat20-29  0.027511   0.102133   0.269  0.78765   
## age_cat30-49  0.063764   0.113771   0.560  0.57517   
## age_cat50-69 -0.387889   0.259240  -1.496  0.13459   
## age_cat70+   -0.639203   0.915770  -0.698  0.48518   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5712.4  on 4166  degrees of freedom
## Residual deviance: 5705.1  on 4159  degrees of freedom
## AIC: 5721.1
## 
## Number of Fisher Scoring iterations: 4

Để thay đổi giá trị tham chiếu của một biến Factor và chuyển giá trị mong muốn lên vị trí đầu tiên, dùng hàm fct_relevel() (xem chương Factors). Ở ví dụ bên dưới, chúng tôi lấy biến age_cat và đặt nhóm tuổi “20-29” làm giá trị tham chiếu trước khi chuyển số liệu đã sửa đổi vào hàm glm().

linelist %>% 
  mutate(age_cat = fct_relevel(age_cat, "20-29", after = 0)) %>% 
  glm(formula = outcome ~ age_cat, family = "binomial") %>% 
  summary()

## 
## Call:
## glm(formula = outcome ~ age_cat, family = "binomial", data = .)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.339  -1.278   1.024   1.080   1.354  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.26125    0.07163   3.647 0.000265 ***
## age_cat0-4   -0.02751    0.10213  -0.269 0.787652    
## age_cat5-9   -0.09041    0.10090  -0.896 0.370220    
## age_cat10-14  0.11069    0.10639   1.040 0.298133    
## age_cat15-19 -0.03308    0.11259  -0.294 0.768934    
## age_cat30-49  0.03625    0.11302   0.321 0.748390    
## age_cat50-69 -0.41540    0.25891  -1.604 0.108625    
## age_cat70+   -0.66671    0.91568  -0.728 0.466546    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5712.4  on 4166  degrees of freedom
## Residual deviance: 5705.1  on 4159  degrees of freedom
## AIC: 5721.1
## 
## Number of Fisher Scoring iterations: 4

In kết quả

Đối với hầu hết các mục đích sử dụng, kết quả đầu ra cần phải có một số sửa đổi. Hàm làm gọn tidy() trong package broom có những tiện lợi để hiển thị kết quả của mô hình.

Ở đây, chúng tôi trình bày cách để kết hợp các kết quả đầu ra của mô hình vào trong một bảng.

1. Lấy lũy thừa logarit của ước lượng tỉ số chênh OR và khoảng tin cậy bằng cách đưa mô hình vào hàm tidy() và thiết lập lũy thừa exponentiate = TRUE và conf.int = TRUE.

model <- glm(outcome ~ age_cat, family = "binomial", data = linelist) %>% 
  tidy(exponentiate = TRUE, conf.int = TRUE) %>%        # exponentiate and produce CIs
  mutate(across(where(is.numeric), round, digits = 2))  # round all numeric columns

Bên dưới là bảng kết quả đầu ra của model:

2. Kết hợp các kết quả của mô hình vào trong một bảng đếm. Dưới đây, chúng tôi tạo một bảng đếm bằng hàm tabyl() từ package janitor, như được đề cập trong chương Bảng mô tả.

counts_table <- linelist %>% 
  janitor::tabyl(age_cat, outcome)

Đây là cách mà bảng counts_table được hiển thị:

Bây giờ chúng ta có thể nối bảng counts_table và kết quả của mô hình model lại với nhau theo chiều ngang bằng hàm nối cột bind_cols() (dplyr). Hãy nhớ rằng đối với hàm bind_cols() các hàng trong hai cấu trúc dữ liệu trên phải được căn chỉnh hoàn hảo. Trong đoạn code này, bởi vì chúng ta đang thực hiện một chuỗi các thuật toán pipe, chúng ta sử dụng dấu . để đại diện cho đối tượng được nối trong bảng đếm counts_table khi chúng tôi nối nó với kết quả mô hình model. Để kết thúc quy trình này, chúng ta sử dụng hàm select() để chọn các cột mong muốn và thứ tự của nó, và cuối cùng áp dụng hàm round() trong base R để làm tròn với hai chữ số thập phân cho tất cả các cột.

combined <- counts_table %>%           # begin with table of counts
  bind_cols(., model) %>%              # combine with the outputs of the regression 
  select(term, 2:3, estimate,          # select and re-order cols
         conf.low, conf.high, p.value) %>% 
  mutate(across(where(is.numeric), round, digits = 2)) ## round to 2 decimal places

Đây là hiển thị của cấu trúc đã được kết hợp, nó được xuất gọn gẽ dưới dạng một hình bằng thông qua một hàm trong package flextable. Chương Trình bày bảng giải thích cách tùy chỉnh các bảng như vậy bằng flextable, hoặc có thể sử dụng các gói lệnh khác như knitr hoặc GT.

combined <- combined %>% 
  flextable::qflextable()

Vòng lặp cho nhiều mô hình đơn biến

Sau đây chúng tôi trình bày một phương pháp sử dụng glm() và tidy() để có một cách tiếp cận đơn giản hơn, xem thêm ở phần gtsummary.

Để thực hiện các mô hình cho một số biến giải thích và cho ra các tỷ số chênh trong phân tích đơn biến (nghĩa là không có kiểm soát lẫn nhau), chúng ta có thể sử dụng các cách tiếp cận dưới đây. Sử dụng hàm str_c() từ package stringr để tạo ra các công thức cho phân tích đơn biến (xem chương Ký tự và chuỗi), thực hiện hàm glm() cho mỗi công thức, chuyển mỗi kết quả đầu ra của glm() đến hàm tidy() và cuối cùng thu gọn lại tất các kết quả đầu ra của mô hình bằng hàm nối dòng bind_rows() từ tidyr. Phương pháp này sử dụng hàm map() từ package purrr để lặp - xem chương [Lặp, vòng lặp và danh sách] để biết thêm thông tin về công cụ này.

1. Tạo một véctơ tên các cột của biến giải thích. Chúng ta đã tạo biến này explanatory_vars trong phần chuẩn bị của chương này.

2. Sử dụng hàm str_c() để tạo các công thức chuỗi với biến kết cục outcome ở bên trái và tên một cột của véctơ explanatory_vars ở bên phải. Dấu chấm . trong hàm này thay thế cho tên cột trong véctơ explanatory_vars.

explanatory_vars %>% str_c("outcome ~ ", .)

## [1] "outcome ~ gender"  "outcome ~ fever"   "outcome ~ chills"  "outcome ~ cough"   "outcome ~ aches"   "outcome ~ vomit"   "outcome ~ age_cat"

3. Đưa các công thức chuỗi này vào hàm map() và đặt ~glm() làm hàm áp dụng cho mỗi đầu vào. Bên trong hàm glm(), thiết lập công thức hồi quy as.formula(.x) trong đó .x sẽ được thay thế bằng các công thức chuỗi đã được tạo bên trên. Hàm map() sẽ lặp từng công thức chuỗi và thực hiện hồi quy cho từng công thức.

4. Kết quả đầu ra của hàm map() đầu tiên sẽ được chuyển đến hàm map() thứ hai mà sử dụng hàm tidy() để làm gọn các kết quả đầu ra.

5. Cuối cùng, kết quả đầu ra của hàm map() thứ hai (một danh sách các data frames đã được làm gọn) được tóm tắt bằng hàm nối dòng bind_rows(), kết quả cho ra một data frame với tất cả các kết quả đơn biến.

models <- explanatory_vars %>%       # begin with variables of interest
  str_c("outcome ~ ", .) %>%         # combine each variable into formula ("outcome ~ variable of interest")
  
  # iterate through each univariate formula
  map(                               
    .f = ~glm(                       # pass the formulas one-by-one to glm()
      formula = as.formula(.x),      # within glm(), the string formula is .x
      family = "binomial",           # specify type of glm (logistic)
      data = linelist)) %>%          # dataset
  
  # tidy up each of the glm regression outputs from above
  map(
    .f = ~tidy(
      .x, 
      exponentiate = TRUE,           # exponentiate 
      conf.int = TRUE)) %>%          # return confidence intervals
  
  # collapse the list of regression outputs in to one data frame
  bind_rows() %>% 
  
  # round all numeric columns
  mutate(across(where(is.numeric), round, digits = 2))

Lúc này, kết quả xuất ra của models dài hơn bởi vì kết quả bây giờ bao gồm các kết quả đầu ra của một số hồi quy đơn biến. Nhấp nút tiếp theo để xem tất cả các hàng của model.

Như lúc trước, chúng ta có thể tạo một bảng đếm từ bộ số liệu linelist cho mỗi biến giải thích, gắn chúng với models, và tạo ra một bảng đẹp. Chúng ta bắt đầu với các biến giải thích này, và lặp lại các biến này thông qua hàm map(). Chúng ta lặp lại qua một hàm do người dùng tạo ra mà liên quan đến việc tạo ra một bảng đếm bằng cách dùng các hàm trong package dplyr Sau đó, kết quả được kết nối trình tự với kết quả của mô hình models.

## for each explanatory variable
univ_tab_base <- explanatory_vars %>% 
  map(.f = 
    ~{linelist %>%                ## begin with linelist
        group_by(outcome) %>%     ## group data set by outcome
        count(.data[[.x]]) %>%    ## produce counts for variable of interest
        pivot_wider(              ## spread to wide format (as in cross-tabulation)
          names_from = outcome,
          values_from = n) %>% 
        drop_na(.data[[.x]]) %>%         ## drop rows with missings
        rename("variable" = .x) %>%      ## change variable of interest column to "variable"
        mutate(variable = as.character(variable))} ## convert to character, else non-dichotomous (categorical) variables come out as factor and cant be merged
      ) %>% 
  
  ## collapse the list of count outputs in to one data frame
  bind_rows() %>% 
  
  ## merge with the outputs of the regression 
  bind_cols(., models) %>% 
  
  ## only keep columns interested in 
  select(term, 2:3, estimate, conf.low, conf.high, p.value) %>% 
  
  ## round decimal places
  mutate(across(where(is.numeric), round, digits = 2))

Bên dưới là cấu trúc số liệu kết nối được tạo ra. Xem chương Trình bày bảng để có thêm ý tưởng về cách chuyển đổi bảng số liệu này thành một bảng đẹp trên HTML (ví dụ như với package flextable).

Xây dựng mô hình

Chúng ta có thể xây dựng mô hình theo từng bước, lưu các mô hình khác nhau với một số biến giải thích. Chúng ta có thể sử dụng kiểm định tỷ số khả dĩ (likelihood-ratio tests) để so sánh các mô hình này bằng cách sử dụng hàm lrtest() từ package lmtest, như dưới đây:

CHÚ Ý: Sử dụng hàn anova(model1, model2, test = "Chisq) trong base R cũng cho kết quả tương tự

model1 <- glm(outcome ~ age_cat, family = "binomial", data = linelist)
model2 <- glm(outcome ~ age_cat + gender, family = "binomial", data = linelist)

lmtest::lrtest(model1, model2)

## Likelihood ratio test
## 
## Model 1: outcome ~ age_cat
## Model 2: outcome ~ age_cat + gender
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   8 -2852.6                     
## 2   9 -2852.6  1 0.0002     0.9883

Một tùy chọn khác là lấy đối tượng của mô hình và sử dụng hàm step() từ package stats. Chỉ rõ hướng lựa chọn biến mà chúng ta muốn sử dụng khi xây dựng mô hình.

## choose a model using forward selection based on AIC
## you can also do "backward" or "both" by adjusting the direction
final_mv_reg <- mv_reg %>%
  step(direction = "forward", trace = FALSE)

Để hiển thị rõ số, chúng ta có thể tắt ký hiệu khoa học trong R bằng lệnh sau

options(scipen=999)

Như được mô tả trong phần phân tích đơn biến, chuyển kết quả đầu ra của mô hình vào hàm tidy() để lấy lũy thừa cho các hệ số và khoảng tin cậy (CIs). Cuối cùng, làm tròn tất cả các cột số với hai số thập phân. Kéo qua để xem tất cả các hàng.

mv_tab_base <- final_mv_reg %>% 
  broom::tidy(exponentiate = TRUE, conf.int = TRUE) %>%  ## get a tidy dataframe of estimates 
  mutate(across(where(is.numeric), round, digits = 2))          ## round 

Đây là hiển thị kết quả dưới dạng data frame looks:

Gộp bằng package gtsummary

Hàm tbl_regression() trong package gtsummary sẽ lấy kết quả đầu tra từ một hồi quy (hàm glm() trong trường hợp này) và tạo ra một bảng tóm tắt đẹp.

## show results table of final regression 
mv_tab <- tbl_regression(final_mv_reg, exponentiate = TRUE)

Hãy xem bảng sau:

mv_tab

Chúng ta cũng có thể kết hợp một số bảng kết quả đầu ra bằng cách dùng hàm tbl_merge() trong package gtsummary. Bây giờ chúng ta hộp các kết quả đa biến với kết quả đơn biến đã được tạo bên trên bằng package gtsummary:

## combine with univariate results 
tbl_merge(
  tbls = list(univ_tab, mv_tab),                          # combine
  tab_spanner = c("**Univariate**", "**Multivariable**")) # set header names

Gộp bằng package dplyr

Một cách khác để gộp các kết quả đơn biến và đa biến từ các hàm glm()/tidy() bằng cách sử dụng các hàm kết nối từ package dplyr.

• Kết nối kết quả đơn biến trước đó (univ_tab_base, chứa được các cột đếm) với kết quả đa biến đã được làm gọn mv_tab_base
  
• Sử dụng hàm select() để giữ lại, sắp xếp lại thứ tự và đặt lại tên các cột mà chúng ta muốn
  
• Sử dụng hàm round() để làm tròn tất cả các cột với hai số thập phân

## combine univariate and multivariable tables 
left_join(univ_tab_base, mv_tab_base, by = "term") %>% 
  ## choose columns and rename them
  select( # new name =  old name
    "characteristic" = term, 
    "recovered"      = "0", 
    "dead"           = "1", 
    "univ_or"        = estimate.x, 
    "univ_ci_low"    = conf.low.x, 
    "univ_ci_high"   = conf.high.x,
    "univ_pval"      = p.value.x, 
    "mv_or"          = estimate.y, 
    "mvv_ci_low"     = conf.low.y, 
    "mv_ci_high"     = conf.high.y,
    "mv_pval"        = p.value.y 
  ) %>% 
  mutate(across(where(is.double), round, 2))   

## # A tibble: 20 x 11
##    characteristic recovered  dead univ_or univ_ci_low univ_ci_high univ_pval mv_or mvv_ci_low mv_ci_high mv_pval
##    <chr>              <dbl> <dbl>   <dbl>       <dbl>        <dbl>     <dbl> <dbl>      <dbl>      <dbl>   <dbl>
##  1 (Intercept)          909  1168    1.28        1.18         1.4       0     1.07       0.83       1.39    0.6 
##  2 gender               916  1174    1           0.88         1.13      0.97  1          0.88       1.14    0.97
##  3 (Intercept)          340   436    1.28        1.11         1.48      0     1.07       0.83       1.39    0.6 
##  4 fever               1485  1906    1           0.85         1.17      0.99  1          0.86       1.18    0.96
##  5 (Intercept)         1472  1877    1.28        1.19         1.37      0     1.07       0.83       1.39    0.6 
##  6 chills               353   465    1.03        0.89         1.21      0.68  1.03       0.89       1.21    0.67
##  7 (Intercept)          272   309    1.14        0.97         1.34      0.13  1.07       0.83       1.39    0.6 
##  8 cough               1553  2033    1.15        0.97         1.37      0.11  1.15       0.96       1.37    0.12
##  9 (Intercept)         1636  2114    1.29        1.21         1.38      0     1.07       0.83       1.39    0.6 
## 10 aches                189   228    0.93        0.76         1.14      0.51  0.93       0.76       1.14    0.5 
## 11 (Intercept)          931  1144    1.23        1.13         1.34      0     1.07       0.83       1.39    0.6 
## 12 vomit                894  1198    1.09        0.96         1.23      0.17  1.09       0.96       1.23    0.17
## 13 (Intercept)          338   427    1.26        1.1          1.46      0     1.07       0.83       1.39    0.6 
## 14 age_cat5-9           365   433    0.94        0.77         1.15      0.54  0.94       0.77       1.15    0.53
## 15 age_cat10-14         273   396    1.15        0.93         1.42      0.2   1.15       0.93       1.41    0.2 
## 16 age_cat15-19         238   299    0.99        0.8          1.24      0.96  0.99       0.79       1.24    0.92
## 17 age_cat20-29         345   448    1.03        0.84         1.26      0.79  1.03       0.84       1.26    0.8 
## 18 age_cat30-49         228   307    1.07        0.85         1.33      0.58  1.06       0.85       1.33    0.61
## 19 age_cat50-69          35    30    0.68        0.41         1.13      0.13  0.68       0.4        1.13    0.14
## 20 age_cat70+             3     2    0.53        0.07         3.2       0.49  0.52       0.07       3.19    0.48