{"id":13527,"date":"2025-09-19T13:27:32","date_gmt":"2025-09-19T13:27:32","guid":{"rendered":"https:\/\/www.mmj.mw\/?p=13527"},"modified":"2025-09-19T13:30:00","modified_gmt":"2025-09-19T13:30:00","slug":"a-practical-guide-to-key-considerations-in-logistic-regression-for-clinical-and-public-health-research-r-tutorial","status":"publish","type":"post","link":"https:\/\/www.mmj.mw\/?p=13527","title":{"rendered":"A Practical Guide to Key Considerations in Logistic Regression for Clinical and Public Health Research: R tutorial"},"content":{"rendered":"\n

Wingston Felix Ng\u2019ambi1<\/sup>, Cosmas Zyambo, Adamson Sinjani Muula3-5<\/sup><\/strong><\/p>\n\n\n\n

  1. Health Economics and Policy Unit, Department of Health Systems and Policy, Kamuzu University of Health Sciences, Lilongwe, Malawi<\/li>
  2. Africa Centre of Excellence in Public Health and Herbal Medicine (ACEPHEM), Kamuzu University of Health Sciences, Blantyre, Malawi<\/li>
  3. Department of Public Health and Family Medicine, University of Zambia, Lusaka, Zambia<\/li>
  4. Professor and Head, Department of Community and Environmental Health, School of Global and Public Health<\/li>
  5. President, ECSA College of Public Health Physicians, Arusha, Tanzania<\/li><\/ol>\n\n\n\n

    Abstract<\/span><\/strong><\/a><\/h4>\n\n\n\n

    Logistic regression is one of the most widely used statistical methods in clinical and public health research, especially for analysing binary outcomes such as disease presence, treatment uptake, or health behaviour change. While its use appears straightforward, important methodological decisions made during data preparation, model specification, and interpretation can substantially influence findings and their policy relevance. This paper offers a practical, step-by-step guide to key considerations in logistic regression, tailored for clinical and public health researchers. Drawing from applied examples in infectious diseases, non-communicable diseases, and health services research, we highlight common pitfalls and best practices for variable selection, model diagnostics, handling missing data, assessing fit, and presenting results. The guide aims to bridge the gap between statistical theory and real-world application, equipping researchers to produce more robust, transparent, and policy-relevant evidence.<\/p>\n\n\n\n

    Key Words<\/span><\/strong><\/a><\/h4>\n\n\n\n

    Logistic regression, Clinical research, Public Health, infectious diseases, non-communicable diseases, health services research<\/p>\n\n\n\n

    Introduction<\/strong><\/h4>\n\n\n\n

    Logistic regression is one of the most widely used statistical tools in clinical and public health research1<\/a><\/sup>. Its versatility, interpretability, and capacity to model binary outcomes make it an indispensable method in epidemiology, health services research, and health policy analysis3<\/a><\/sup>. In clinical and public health, binary outcomes are common8<\/a><\/sup>; for example, whether an individual is vaccinated, whether a patient adheres to treatment, or whether a pregnant woman delivers in a health facility. Logistic regression allows researchers to examine the relationship between such outcomes and multiple predictors, while adjusting for potential confounding factors9<\/a><\/sup>. In global health research contexts, logistic regression underpins analyses from flagship datasets such as the Demographic and Health Surveys (DHS), the Multiple Indicator Cluster Surveys (MICS), the World Health Organization\u2019s STEPS surveys, and disease-specific surveys such as the Population-based HIV Impact Assessments (PHIA). These data sources are often the backbone of national health policy reviews, donor-funded programme evaluations, and academic studies 10<\/a><\/sup>. For instance, logistic regression has been used to identify determinants of HIV testing uptake across sub-Saharan Africa, factors influencing childhood immunisation in South Asia, and socio-economic predictors of maternal health service use in Latin America. Despite its broad adoption, the method is often misapplied or under-reported. Studies may omit critical details about variable selection, ignore checks for key assumptions such as linearity in the logit, or misinterpret odds ratios as risk ratios. Where many datasets are collected using complex survey designs, a frequent oversight is failing to account for clustering, stratification, and sampling weights11<\/a><\/sup>. This can lead to biased estimates, misleading standard errors, and ultimately, practice and policy recommendations that are not supported by the data12<\/a><\/sup>.<\/p>\n\n\n\n

    The gap between statistical theory and its practical application is a constant threat in understanding research findings. Researchers may be constrained by small sample sizes13<\/a><\/sup>, incomplete data14<\/a><\/sup>, or limited access to statistical training and software15<\/a><\/sup>. These challenges can result in oversimplified analyses that miss important interactions or produce unstable estimates16<\/a><\/sup>. For example, a study examining determinants of hypertension treatment in a country may exclude socio-economic variables because of missing data, or might not test for interaction effects between age and sex, thereby overlooking potentially important policy insights. The aim of this paper is to provide a practical, step-by-step guide to key considerations when applying the logistic regression in public health research, with a specific focus on contexts relevant to both global and African health studies. This guide bridges the gap between statistical best practice and the realities of working with large-scale survey data, health facility datasets, and programmatic monitoring systems. It is written for applied researchers, programme evaluators, postgraduate students, and public health practitioners who may not be statisticians by training but who regularly use logistic regression to generate evidence for decision-making.<\/p>\n\n\n\n

    Data sets analysed<\/strong><\/p>\n\n\n\n

    Throughout the paper, we combine statistical principles with worked examples in R, using real-world datasets from both global and African sources. Examples include:<\/p>\n\n\n\n

    – Predictors of full immunisation among children aged 12\u201323 months using the Malawi Demographic and Health Survey (DHS) 2015\u201316 data MDHS 2015-2016<\/a>. Below is the code that creates the dhs_mw dataset for analysis.<\/p>\n\n\n\n

    # Always set working directory to the current project folder
    setwd(getwd())

    # Load required packages ———————————————

    library(haven)      # for reading Stata (.dta) files like DHS datasets
    library(dplyr)      # for data wrangling (filtering, mutating, selecting)<\/p>\n\n\n\n


    Attaching package: ‘dplyr’<\/p>\n\n\n\n

    The following objects are masked from ‘package:stats’:

        filter, lag<\/p>\n\n\n\n

    The following objects are masked from ‘package:base’:

        intersect, setdiff, setequal, union<\/p>\n\n\n\n

    library(labelled)   # for handling variable\/value labels from DHS
    library(mice)       # for multiple imputation of missing values<\/p>\n\n\n\n


    Attaching package: ‘mice’<\/p>\n\n\n\n

    The following object is masked from ‘package:stats’:

        filter<\/p>\n\n\n\n

    The following objects are masked from ‘package:base’:

        cbind, rbind<\/p>\n\n\n\n

    library(survey)     # for survey-weighted analyses (logistic regression)<\/p>\n\n\n\n

    Loading required package: grid<\/p>\n\n\n\n

    Loading required package: Matrix<\/p>\n\n\n\n

    Loading required package: survival<\/p>\n\n\n\n


    Attaching package: ‘survey’<\/p>\n\n\n\n

    The following object is masked from ‘package:graphics’:

        dotchart<\/p>\n\n\n\n

    # 1. Load DHS dataset ————————————————

    dhs_mw <- read_dta(“MWKR7AFL.DTA”)
    # Reads the DHS Malawi Kids Recode file (KR) in Stata format (.dta).
    # The object ‘dhs_mw’ now contains the full DHS dataset for analysis.

    # 2. Select children age 12\u201323 months ——————————–

    dhs_mw <- dhs_mw %>%
    \u00a0 filter(b8 >= 1 & b8 <= 2)
    # b8 = age of child in years
    # We restrict the sample to children 12\u201323 months old, since they are
    # expected to have completed all basic vaccinations.

    # 3. Define full immunisation outcome ——————————–

    dhs_mw <- dhs_mw %>%
    \u00a0 mutate(
    \u00a0\u00a0\u00a0 full_immunisation = case_when(
    \u00a0\u00a0\u00a0\u00a0\u00a0 # child fully immunised if all required vaccines are 1,2,3
    \u00a0\u00a0\u00a0\u00a0\u00a0 h2 %in% c(1,2,3) &\u00a0\u00a0 # BCG
    \u00a0\u00a0\u00a0\u00a0\u00a0 h3 %in% c(1,2,3) &\u00a0\u00a0 # DPT1
    \u00a0\u00a0\u00a0\u00a0\u00a0 h4 %in% c(1,2,3) &\u00a0\u00a0 # DPT2
    \u00a0\u00a0\u00a0\u00a0\u00a0 h5 %in% c(1,2,3) &\u00a0\u00a0 # DPT3
    \u00a0\u00a0\u00a0\u00a0\u00a0 h6 %in% c(1,2,3) &\u00a0\u00a0 # Polio1
    \u00a0\u00a0\u00a0\u00a0\u00a0 h7 %in% c(1,2,3) &\u00a0\u00a0 # Polio2
    \u00a0\u00a0\u00a0\u00a0\u00a0 h8 %in% c(1,2,3) &\u00a0\u00a0 # Polio3
    \u00a0\u00a0\u00a0\u00a0\u00a0 h9 %in% c(1,2,3) &\u00a0\u00a0 # Measles
    \u00a0\u00a0\u00a0\u00a0\u00a0 h9a %in% c(1,2,3)\u00a0\u00a0\u00a0 # Sometimes measles2 \/ rubella etc depending on DHS round
    \u00a0\u00a0\u00a0\u00a0\u00a0 ~ 1,
    \u00a0\u00a0\u00a0\u00a0\u00a0 TRUE ~ 0
    \u00a0\u00a0\u00a0 )
    \u00a0 )
    # Creates a binary outcome variable ‘full_immunisation’.
    # A child is coded 1 if they received all required vaccines:
    #\u00a0\u00a0 BCG, 3 doses of DPT, 3 doses of polio, and measles.
    # Otherwise, coded as 0.
    # NOTE: variable names (bcg, dpt1, etc.) may differ in your DHS file
    # and need to be checked.

    # 4. Recode predictors and keep analysis variables ——————-

    dhs_mw <- dhs_mw %>%
    \u00a0 mutate(
    \u00a0\u00a0\u00a0 child_age = b8,\u00a0\u00a0 # keep child age in months as continuous predictor
    \u00a0\u00a0\u00a0
    \u00a0\u00a0\u00a0 mother_education = factor(v106,\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 # mother’s education level
    \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 levels = c(0,1,2,3),
    \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 labels = c(“none”,”primary”,”secondary”,”higher”)),
    \u00a0\u00a0\u00a0
    \u00a0\u00a0\u00a0 wealth_quintile = factor(v190,\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 # household wealth quintile
    \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 levels = 1:5,
    \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 labels = c(“poorest”,”poorer”,”middle”,”richer”,”richest”)),
    \u00a0\u00a0\u00a0
    \u00a0\u00a0 \u00a0urban = ifelse(v025 == 1, 1, 0),\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 # place of residence
    \u00a0\u00a0\u00a0 # v025 = type of residence (1=urban, 2=rural).
    \u00a0\u00a0\u00a0 # We recode to binary (1=urban, 0=rural).
    \u00a0\u00a0\u00a0
    \u00a0\u00a0\u00a0 mother_age = factor(v013,\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 # mother’s age group
    \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 levels = 1:7,
    \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 labels = c(“15-19″,”20-24″,”25-29”,
    \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 “30-34″,”35-39″,”40-44″,”45-49”)),
    \u00a0\u00a0\u00a0
    \u00a0\u00a0\u00a0 region = factor(v024,\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 # region of residence
    \u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0levels = c(1,2,3),
    \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 labels = c(“Northern”,”Central”,”Southern”)),
    \u00a0\u00a0\u00a0
    \u00a0\u00a0\u00a0 cluster = v021,\u00a0\u00a0\u00a0 # cluster (primary sampling unit, PSU)
    \u00a0\u00a0\u00a0 strata\u00a0 = v022,\u00a0\u00a0\u00a0 # survey strata
    \u00a0\u00a0\u00a0 weight\u00a0 = v005 \/ 1000000\u00a0\u00a0 # sample weight (normalize by 1,000,000 as per DHS)
    \u00a0 ) %>%
    \u00a0 select(full_immunisation, child_age, mother_education, mother_age,
    \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 wealth_quintile, urban, region, cluster, strata, weight)
    # After recoding, we keep only the variables relevant for the analysis.

    #Dealing with data types so that we are able to run all the commands perfectly
    dhs_mw$mother_education <- factor(dhs_mw$mother_education,
    \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 levels = c(“none”, “primary”, “secondary”, “higher”))
    dhs_mw$child_age <- factor(dhs_mw$child_age, levels = c(1,2))
    dhs_mw$urban <- factor(dhs_mw$urban, levels = c(0,1))
    dhs_mw$wealth_quintile <- factor(dhs_mw$wealth_quintile,
    \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 levels = c(“poorest”, “poorer”, “middle”, “richer”, “richest”))<\/p>\n\n\n\n

    By grounding the discussion in real datasets, we demonstrate how key considerations such as outcome coding, variable selection, model diagnostics, and interpretation translate into practical steps. We also highlight common pitfalls and provide recommendations to ensure logistic regression results are robust, transparent, and practice and policy-relevant. In doing so, this paper aims to demystify logistic regression for applied public health research, contexts where data complexity and policy urgency demand both statistical rigour and pragmatic decision-making.<\/p>\n\n\n\n

    Step-by-step guide to key considerations<\/strong><\/p>\n\n\n\n

    This section gives a hands-on, practical walk through the most important choices you will make when using logistic regression in clinical and public health work. Each part explains why the step matters, shows common mistakes, and gives short R code you can copy and adapt. We use short, clear language, one idea per sentence, and real public health examples from Africa and elsewhere.<\/p>\n\n\n\n

    Step 1: Define the research question<\/em> A clear question is the most important start. Ask, what is the exact outcome. Ask, is it naturally binary, or am I forcing a binary split. Be explicit about whether your aim is explanation, causal inference, or prediction17<\/a><\/sup>. Each aim needs a slightly different approach. A precise research question is the foundation of any analysis. This requires clear specification of the outcome, with attention to whether it is naturally binary or whether dichotomisation is being imposed on a continuous or ordinal variable, which can reduce information18<\/a><\/sup>. Equally important is explicit definition of the analytic aim. Explanatory analyses prioritise identifying associations and patterns between variables, whereas causal inference focuses on estimating the effect of an exposure or intervention under defined assumptions19<\/a><\/sup>. In contrast, prediction models aim to maximise accuracy in forecasting outcomes, often using machine learning methods that prioritise performance over interpretability. Distinguishing between these aims at the outset ensures alignment between research design, choice of methods, and interpretation of results20<\/a><\/sup>.<\/p>\n\n\n\n

    Example, explanatory question<\/em>: \u201cWhat factors predict facility delivery in rural Malawi?\u201d Here you want effect estimates that adjust for confounding. Use odds ratios carefully, and consider risk differences for policy clarity.<\/p>\n\n\n\n

    Example, predictive question<\/em>: \u201cCan we predict which children will be fully immunised in Malawi?\u201d Here you focus on discrimination and calibration, you care less about causal interpretation.<\/p>\n\n\n\n

    Practical checks, before modelling<\/em>: Before starting any modelling, it is important to carry out some practical checks to ensure your analysis is robust. First, inspect the frequency of your outcome: very rare events or extremely common outcomes can pose challenges for standard regression methods and may require alternative approaches or careful interpretation. Second, decide on your reference category for categorical variables, code it clearly (for example, 0 versus 1), and document this choice explicitly in your methods section. Doing so improves transparency, reproducibility, and clarity for readers, reviewers, and policymakers who may rely on your findings.<\/p>\n\n\n\n

    R, quick checks:<\/p>\n\n\n\n

    # check outcome distribution
    table(dhs_mw$full_immunisation)<\/p>\n\n\n\n


       0    1
    5626  874<\/p>\n\n\n\n

    prop.table(table(dhs_mw$full_immunisation))<\/p>\n\n\n\n


            0         1
    0.8655385 0.1344615<\/p>\n\n\n\n

    Step 2: Data preparation<\/em> Outcome coding<\/em>: This section addresses outcome coding, predictor coding, handling of missing data, survey design features, and simple data transformations. For outcome coding, clarity is essential. Outcomes should be coded consistently, often using 0 and 1 to indicate absence and presence of the event, respectively3<\/a><\/sup>. It is critical to state explicitly which value represents the event of interest. For example, in the Malawi Demographic and Health Survey (DHS) child immunisation dataset, the outcome can be coded as 1 = fully immunised and 0 = not fully immunised. Such explicit coding not only reduces ambiguity but also facilitates reproducibility and comparability across analyses. R, example:<\/p>\n\n\n\n

    table(dhs_mw$full_immunisation, useNA = “ifany”)<\/p>\n\n\n\n


       0    1
    5626  874<\/p>\n\n\n\n

    Predictor coding<\/em>: For predictor coding, simplicity and interpretability should guide decisions. Categories should be meaningful and analytically relevant, while avoiding excessive fragmentation into many small groups with very few observations, which can compromise statistical power and model stability21<\/a><\/sup>. When continuous variables are used in interaction terms, centring them around their mean is recommended, as this improves the interpretability of regression coefficients and reduces multicollinearity22<\/a><\/sup>.<\/p>\n\n\n\n

    Missing data<\/em>: For missing data, exploration of patterns is the first step. It is important to assess whether data are missing completely at random (MCAR), missing at random (MAR), or missing not at random (MNAR), as each mechanism has different implications for analysis14<\/a><\/sup>. In many health survey settings, multiple imputation provides a robust approach for handling missingness and yields less biased estimates compared to ad hoc methods25<\/a><\/sup>. Complete-case analysis should be avoided when the proportion of missing data is greater than 10%, as this can reduce statistical power and lead to biased results if the missingness is not completely at random26<\/a><\/sup>.<\/p>\n\n\n\n

    R, example with mice:<\/p>\n\n\n\n

    library(mice)
    # md.pattern(dhs_mw) # Our data is completely observed with no missing data
    imp <- mice(dhs_mw %>% select(full_immunisation, child_age, mother_education, wealth_quintile, urban),
                m = 5, method = ‘pmm’, seed = 123)<\/p>\n\n\n\n


     iter imp variable
      1   1
      1   2
      1   3
      1   4
      1   5
      2   1
      2   2
      2   3
      2   4
      2   5
      3   1
      3   2
      3   3
      3   4
      3   5
      4   1
      4   2
      4   3
      4   4
      4   5
      5   1
      5   2
      5   3
      5   4
      5   5<\/p>\n\n\n\n

    mod_imp <- with(imp, glm(full_immunisation ~ child_age + mother_education + wealth_quintile + urban,
                             family = binomial))
    pool(mod_imp)<\/p>\n\n\n\n

    Class: mipo    m = 5
                            term m    estimate        ubar b           t dfcom
    1                (Intercept) 5 -2.02286323 0.017112417 0 0.017112417  6490
    2                 child_age2 5 -0.17364443 0.005336868 0 0.005336868  6490
    3    mother_educationprimary 5  0.04161282 0.014670901 0 0.014670901  6490
    4  mother_educationsecondary 5  0.11769675 0.021149348 0 0.021149348  6490
    5     mother_educationhigher 5 -0.03311987 0.094604028 0 0.094604028  6490
    6      wealth_quintilepoorer 5  0.13756188 0.013136974 0 0.013136974  6490
    7      wealth_quintilemiddle 5  0.17057447 0.013670904 0 0.013670904  6490
    8      wealth_quintilericher 5  0.31267124 0.014138052 0 0.014138052  6490
    9     wealth_quintilerichest 5  0.29185892 0.019885581 0 0.019885581  6490
    10                    urban1 5  0.10163447 0.012979206 0 0.012979206  6490
             df riv lambda          fmi
    1  6487.247   0      0 0.0003081547
    2  6487.247   0      0 0.0003081547
    3  6487.247   0      0 0.0003081547
    4  6487.247   0      0 0.0003081547
    5  6487.247   0      0 0.0003081547
    6  6487.247   0      0 0.0003081547
    7  6487.247   0      0 0.0003081547
    8  6487.247   0      0 0.0003081547
    9  6487.247   0      0 0.0003081547
    10 6487.247   0      0 0.0003081547<\/p>\n\n\n\n

    Survey design and weights<\/em>: Survey design and weights require special attention, as many  national health datasets are based on complex sampling strategies that involve stratification, multi-stage sampling, clustering, and unequal probabilities of selection29<\/a><\/sup>. Ignoring these design features can result in biased standard errors and, in some cases, biased point estimates30<\/a><\/sup>. Analyses should therefore apply survey-aware methods that incorporate weights and account for clustering to ensure valid statistical inference30<\/a><\/sup>.<\/p>\n\n\n\n

    R, using survey package:<\/p>\n\n\n\n

    library(survey)

    dhs_design <- svydesign(ids = ~cluster, strata = ~strata, weights = ~weight, data = dhs_mw, nest = TRUE)

    svyglm(full_immunisation ~ child_age + mother_education + wealth_quintile + urban,
           design = dhs_design, family = quasibinomial())<\/p>\n\n\n\n

    Stratified 1 – level Cluster Sampling design (with replacement)
    With (849) clusters.
    svydesign(ids = ~cluster, strata = ~strata, weights = ~weight,
        data = dhs_mw, nest = TRUE)

    Call:  svyglm(formula = full_immunisation ~ child_age + mother_education +
        wealth_quintile + urban, design = dhs_design, family = quasibinomial())

    Coefficients:
                  (Intercept)                 child_age2 
                     -2.24446                   -0.15802 
      mother_educationprimary  mother_educationsecondary 
                      0.19456                    0.29615 
       mother_educationhigher      wealth_quintilepoorer 
                      0.52557                    0.27186 
        wealth_quintilemiddle      wealth_quintilericher 
                      0.22911                    0.47902 
       wealth_quintilerichest                     urban1 
                      0.33950                    0.05492 

    Degrees of Freedom: 6499 Total (i.e. Null);  784 Residual
    Null Deviance:      5135
    Residual Deviance: 5097     AIC: NA<\/p>\n\n\n\n

    Transformations and outliers<\/em>: For transformations and outliers, distributions of continuous predictors should be carefully examined. Skewed predictors may benefit from log transformations or flexible spline functions to improve model fit and interpretability31<\/a><\/sup>. Extreme values should be handled cautiously; winsorising or truncating should only be applied when there is clear justification, as inappropriate handling of outliers can distort results and reduce external validity32<\/a><\/sup>.<\/p>\n\n\n\n

    R, centring example:<\/p>\n\n\n\n

    dhs_mw <- dhs_mw %>%
      mutate(child_age = case_when(
        child_age == 1 ~ 0,
        child_age == 2 ~ 1,
        TRUE ~ NA_real_  # keep missing or unexpected values as NA
      ))

    # Check the result
    table(dhs_mw$child_age, useNA = “ifany”)<\/p>\n\n\n\n


       0    1
    3248 3252<\/p>\n\n\n\n

    dhs_mw <- dhs_mw %>% mutate(child_age_c = child_age – mean(child_age, na.rm = TRUE))<\/p>\n\n\n\n

    Document every recode<\/em>: Finally, every recoding step should be documented in detail. Transparent documentation supports reproducibility, allows other researchers to replicate the analysis, and helps informed policymakers and reviewers trust the validity of the results33<\/a><\/sup>.<\/p>\n\n\n\n

    Step 3: Variable selection<\/em> Choose predictors with purpose<\/em>: Selection of predictors should be grounded in theory, prior evidence, and where possible, formal tools such as directed acyclic graphs (DAGs) to make assumptions explicit34<\/a><\/sup>. This approach ensures that models reflect plausible causal structures and avoid spurious associations35<\/a><\/sup>. While data-driven methods, such as automated variable selection, can be valuable in prediction tasks, they are less appropriate for explanatory or causal analyses, where reliance on prior knowledge provides stronger justification for including or excluding variables36<\/a><\/sup>. A clear rationale for each predictor strengthens the credibility and interpretability of the final model37<\/a><\/sup>.<\/p>\n\n\n\n

    Events per variable (EPV)<\/em>: The number of outcome events relative to the number of predictors is a key consideration in model building. A widely cited rule of thumb is to maintain at least ten outcome events for each parameter estimated, often referred to as the events-per-variable (EPV) guideline38<\/a><\/sup>. When EPV is low, models risk instability, overfitting, and inflated standard errors. In such situations, strategies include reducing the number of predictors, collapsing sparse categories into broader groups, or applying penalised regression methods such as ridge, lasso, or elastic net. These approaches help preserve model stability while still capturing relevant information39<\/a><\/sup>.<\/p>\n\n\n\n

    R, check EPV:<\/p>\n\n\n\n

    events <- sum(dhs_mw$full_immunisation == 1, na.rm = TRUE)
    num_params <- 6  # approximate number of model coefficients
    epv <- events \/ num_params
    epv<\/p>\n\n\n\n

    1 145.6667<\/p>\n\n\n\n

    Automated stepwise selection is common<\/em>: Automated stepwise selection is a common approach for reducing the number of predictors in a model, but it is inherently unstable and can produce results that are sensitive to small changes in the data. When this method is employed, it is essential to report its use transparently, including the criteria for adding or removing variables8<\/a><\/sup>. To ensure that the final model is robust and reliable, validation techniques such as bootstrapping or cross-validation should be applied, providing an assessment of model stability and predictive performance beyond the original sample40<\/a><\/sup>.<\/p>\n\n\n\n

    R, stepwise warning example:<\/p>\n\n\n\n

    library(MASS)<\/p>\n\n\n\n


    Attaching package: ‘MASS’<\/p>\n\n\n\n

    The following object is masked from ‘package:dplyr’:

        select<\/p>\n\n\n\n

    # List all variables in `dhs_mw`
    for (var in names(dhs_mw)) {
      print(var)
    }<\/p>\n\n\n\n

    1 “full_immunisation”
    1 “child_age”
    1 “mother_education”
    1 “mother_age”
    1 “wealth_quintile”
    1 “urban”
    1 “region”
    1 “cluster”
    1 “strata”
    1 “weight”
    1 “child_age_c”<\/p>\n\n\n\n

    full_mod <- glm(full_immunisation ~., family = binomial, data = dhs_mw)
    step_mod <- stepAIC(full_mod, direction = “both”, trace = FALSE)
    summary(step_mod)<\/p>\n\n\n\n


    Call:
    glm(formula = full_immunisation ~ child_age + wealth_quintile +
        region, family = binomial, data = dhs_mw)

    Coefficients:
                           Estimate Std. Error z value Pr(>|z|)   
    (Intercept)            -2.06314    0.12461 -16.557  < 2e-16 ***
    child_age              -0.16720    0.07325  -2.283 0.022459 * 
    wealth_quintilepoorer   0.14488    0.11461   1.264 0.206198   
    wealth_quintilemiddle   0.21043    0.11700   1.798 0.072102 . 
    wealth_quintilericher   0.36978    0.11733   3.152 0.001624 **
    wealth_quintilerichest  0.40253    0.11774   3.419 0.000629 ***
    regionCentral           0.35952    0.10316   3.485 0.000492 ***
    regionSouthern         -0.17059    0.10419  -1.637 0.101574   
    —
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1

    (Dispersion parameter for binomial family taken to be 1)

        Null deviance: 5132.1  on 6499  degrees of freedom
    Residual deviance: 5068.9  on 6492  degrees of freedom
    AIC: 5084.9

    Number of Fisher Scoring iterations: 4<\/p>\n\n\n\n

    library(broom)
    library(knitr)

    # Assume step_mod is your fitted glm\/logistic model
    # step_mod <- glm(full_immunisation ~ child_age + sex + …, family = binomial(), data = dhs_mw)

    # ——————————-
    # Create a tidy table
    # ——————————-
    tidy_table <- broom::tidy(step_mod) %>%
      mutate(
        OR = exp(estimate),                                  # convert logOR to OR
        lower_CI = exp(estimate – 1.96 * std.error),        # 95% CI lower
        upper_CI = exp(estimate + 1.96 * std.error),        # 95% CI upper
        p_value = p.value                                   # p-value
      ) %>%
      dplyr::select(term, OR, lower_CI, upper_CI, p_value)

    # ——————————-
    # Print nicely
    # ——————————-
    knitr::kable(tidy_table, digits = 3, caption = “Logistic Regression Results: OR and 95% CI”)<\/p>\n\n\n\n

    Logistic Regression Results: OR and 95% CI<\/p>\n\n\n\n

    term<\/td>OR<\/td>lower_CI<\/td>upper_CI<\/td>p_value<\/td><\/tr><\/thead>
    (Intercept)<\/td>0.127<\/td>0.100<\/td>0.162<\/td>0.000<\/td><\/tr>
    child_age<\/td>0.846<\/td>0.733<\/td>0.977<\/td>0.022<\/td><\/tr>
    wealth_quintilepoorer<\/td>1.156<\/td>0.923<\/td>1.447<\/td>0.206<\/td><\/tr>
    wealth_quintilemiddle<\/td>1.234<\/td>0.981<\/td>1.552<\/td>0.072<\/td><\/tr>
    wealth_quintilericher<\/td>1.447<\/td>1.150<\/td>1.822<\/td>0.002<\/td><\/tr>
    wealth_quintilerichest<\/td>1.496<\/td>1.187<\/td>1.884<\/td>0.001<\/td><\/tr>
    regionCentral<\/td>1.433<\/td>1.170<\/td>1.754<\/td>0.000<\/td><\/tr>
    regionSouthern<\/td>0.843<\/td>0.687<\/td>1.034<\/td>0.102<\/td><\/tr><\/tbody><\/table><\/figure>\n\n\n\n

    Penalised regression for many predictors<\/em>: Penalised regression methods are useful when models include many predictors or when the primary goal is accurate prediction39<\/a><\/sup>. Techniques such as lasso regression are widely used because they simultaneously perform variable selection and shrinkage, reducing overfitting and improving model generalisability. By constraining or penalising coefficient estimates, these methods stabilise models with high-dimensional data while retaining the most informative predictors, making them particularly valuable in settings with limited events per variable or complex predictor structures40<\/a><\/sup>.<\/p>\n\n\n\n

    R, lasso example with glmnet:<\/p>\n\n\n\n

    library(glmnet)<\/p>\n\n\n\n

    Loaded glmnet 4.1-8<\/p>\n\n\n\n

    x <- model.matrix(full_immunisation ~ child_age + mother_education + wealth_quintile + urban + region+strata+cluster+mother_age, data = dhs_mw),-1
    y <- dhs_mw$full_immunisation
    cvfit <- cv.glmnet(x, y, family = “binomial”, alpha = 1)
    coef(cvfit, s = “lambda.min”)<\/p>\n\n\n\n

    20 x 1 sparse Matrix of class “dgCMatrix”
                                        s1
    (Intercept)               -1.919637609
    child_age                 -0.104732144
    mother_educationprimary    .         
    mother_educationsecondary  0.050703849
    mother_educationhigher     .         
    wealth_quintilepoorer      .         
    wealth_quintilemiddle      .         
    wealth_quintilericher      0.137566617
    wealth_quintilerichest     0.135410726
    urban1                     0.052631921
    regionCentral              0.287798303
    regionSouthern            -0.166838474
    strata                     .         
    cluster                    .         
    mother_age20-24           -0.002730957
    mother_age25-29            .         
    mother_age30-34            0.031755025
    mother_age35-39            .         
    mother_age40-44            .         
    mother_age45-49           -0.166722734<\/p>\n\n\n\n

    Collinearity<\/em>: Collinearity among predictors can compromise model stability and interpretability8<\/a><\/sup>. It is important to assess correlations using measures such as variance inflation factors (VIF) and to address highly correlated variables by either dropping redundant predictors, combining them into meaningful composite measures, or applying dimensionality reduction techniques such as principal component analysis8<\/a><\/sup>. Managing collinearity ensures that coefficient estimates remain reliable and that the model provides clear and interpretable insights41<\/a><\/sup>.<\/p>\n\n\n\n

    R, VIF:<\/p>\n\n\n\n

    library(car)<\/p>\n\n\n\n

    Loading required package: carData<\/p>\n\n\n\n


    Attaching package: ‘car’<\/p>\n\n\n\n

    The following object is masked from ‘package:dplyr’:

        recode<\/p>\n\n\n\n

    vif(glm(full_immunisation ~ child_age + mother_education + wealth_quintile + urban, family = binomial, data = dhs_mw))<\/p>\n\n\n\n

                         GVIF Df GVIF^(1\/(2*Df))
    child_age        1.002348  1        1.001173
    mother_education 1.315415  3        1.046752
    wealth_quintile  1.698182  4        1.068435
    urban            1.447586  1        1.203157<\/p>\n\n\n\n

    Step 4: Model specification<\/em> Model specification involves selecting the appropriate functional form for predictors, testing potential interactions, and deciding whether to include random effects when the data exhibit clustering beyond the survey design42<\/a><\/sup>. Proper specification ensures that the model accurately captures the underlying relationships without introducing bias or misspecification. Considering interactions allows for more nuanced understanding of how predictors jointly influence the outcome, while random effects account for unobserved heterogeneity in clustered data, improving both the validity and interpretability of the model42<\/a><\/sup>.<\/p>\n\n\n\n

    Linearity in the logit<\/em>: Logistic regression assumes that the log odds of the outcome change linearly with each continuous predictor42<\/a><\/sup>. It is important to assess this assumption, either graphically or using formal statistical tests44<\/a><\/sup>. If the relationship is found to be non-linear, flexible approaches such as spline functions or fractional polynomials can be applied to capture the true shape of the association, ensuring that the model accurately represents the data and improves predictive performance46<\/a><\/sup>.<\/p>\n\n\n\n

    R, simple graphical check and natural spline:<\/p>\n\n\n\n

    library(splines)
    mod_linear <- glm(full_immunisation ~ child_age + mother_education + wealth_quintile + urban, family = binomial, data = dhs_mw)
    mod_spline <- glm(full_immunisation ~ ns(child_age, df = 4) + mother_education + wealth_quintile + urban, family = binomial, data = dhs_mw)<\/p>\n\n\n\n

    Warning in ns(child_age, df = 4): shoving ‘interior’ knots matching boundary
    knots to inside<\/p>\n\n\n\n

    anova(mod_linear, mod_spline, test = “Chisq”)<\/p>\n\n\n\n

    Analysis of Deviance Table

    Model 1: full_immunisation ~ child_age + mother_education + wealth_quintile +
        urban
    Model 2: full_immunisation ~ ns(child_age, df = 4) + mother_education +
        wealth_quintile + urban
      Resid. Df Resid. Dev Df Deviance Pr(>Chi)
    1      6490     5110.5                    
    2      6490     5110.5  0        0        <\/p>\n\n\n\n

    For plotting the logit:<\/p>\n\n\n\n

    dhs_mw <- dhs_mw %>% mutate(pred = predict(mod_linear, type = “response”),
                                logit = log(pred \/ (1 – pred)))
    library(ggplot2)
    ggplot(dhs_mw, aes(x = child_age, y = logit)) + geom_point(alpha = 0.3) + geom_smooth()<\/p>\n\n\n\n

    `geom_smooth()` using method = ‘gam’ and formula = ‘y ~ s(x, bs = “cs”)’<\/p>\n\n\n\n

    Warning: Failed to fit group -1.
    Caused by error in `smooth.construct.cr.smooth.spec()`:
    ! x has insufficient unique values to support 10 knots: reduce k.<\/p>\n\n\n\n

    \"\"<\/figure>\n\n\n\n

    Figure 1: Distribution of full immunization by age of the child<\/strong><\/p>\n\n\n\n

    Interactions<\/em>: Interactions should be included in a model only when there is a plausible rationale based on theory or prior evidence8<\/a><\/sup>. Once included, interactions should be tested statistically and visualised to aid interpretation. In many policy-relevant analyses, certain interactions can be particularly meaningful47<\/a><\/sup>; for example, the effect of education on health outcomes may differ between urban and rural residents. Careful consideration of interactions enhances the model\u2019s ability to capture nuanced relationships and provides more actionable insights for decision-making.<\/p>\n\n\n\n

    R, example interaction:<\/p>\n\n\n\n

    mod_int <- glm(full_immunisation ~ child_age + mother_education * urban + wealth_quintile, family = binomial, data = dhs_mw)
    summary(mod_int)<\/p>\n\n\n\n


    Call:
    glm(formula = full_immunisation ~ child_age + mother_education *
        urban + wealth_quintile, family = binomial, data = dhs_mw)

    Coefficients:
                                      Estimate Std. Error z value Pr(>|z|)   
    (Intercept)                      -1.996732   0.132556 -15.063  < 2e-16 ***
    child_age                        -0.174629   0.073106  -2.389  0.01691 * 
    mother_educationprimary           0.008395   0.124631   0.067  0.94630   
    mother_educationsecondary         0.092991   0.155946   0.596  0.55097   
    mother_educationhigher            0.009672   0.506280   0.019  0.98476   
    urban1                           -0.398505   0.541564  -0.736  0.46183   
    wealth_quintilepoorer             0.140799   0.114704   1.227  0.21964   
    wealth_quintilemiddle             0.172510   0.117111   1.473  0.14074   
    wealth_quintilericher             0.314012   0.119719   2.623  0.00872 **
    wealth_quintilerichest            0.291767   0.141237   2.066  0.03885 * 
    mother_educationprimary:urban1    0.547274   0.555365   0.985  0.32441   
    mother_educationsecondary:urban1  0.495062   0.564461   0.877  0.38046   
    mother_educationhigher:urban1     0.401599   0.796192   0.504  0.61398   
    —
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1

    (Dispersion parameter for binomial family taken to be 1)

        Null deviance: 5132.1  on 6499  degrees of freedom
    Residual deviance: 5109.4  on 6487  degrees of freedom
    AIC: 5135.4

    Number of Fisher Scoring iterations: 4<\/p>\n\n\n\n

    Random effects<\/em>: When data have a hierarchical or clustered structure that is not fully accounted for by survey weights, multilevel models with random effects should be considered48<\/a><\/sup>. For example, children may be nested within clusters, and clusters within districts. Multilevel models explicitly capture between-cluster variation and can provide either cluster-specific or population-averaged estimates, depending on the choice of link function and the interpretation required49<\/a><\/sup>. Incorporating random effects improves model accuracy, accounts for unobserved heterogeneity, and ensures valid inference in clustered data.<\/p>\n\n\n\n

    R, simple multilevel example with lme4:<\/p>\n\n\n\n

    library(lme4)<\/p>\n\n\n\n

    Warning: package ‘lme4’ was built under R version 4.5.1<\/p>\n\n\n\n

    mod_mixed <- glmer(full_immunisation ~ child_age + mother_education + (1 | cluster), family = binomial, data = dhs_mw)
    summary(mod_mixed)<\/p>\n\n\n\n

    Generalized linear mixed model fit by maximum likelihood (Laplace
      Approximation) glmerMod
     Family: binomial  ( logit )
    Formula: full_immunisation ~ child_age + mother_education + (1 | cluster)
       Data: dhs_mw

          AIC       BIC    logLik -2*log(L)  df.resid
       4939.8    4980.4   -2463.9    4927.8      6494

    Scaled residuals:
        Min      1Q  Median      3Q     Max
    -1.0933 -0.3682 -0.2737 -0.2371  3.9218

    Random effects:
     Groups  Name        Variance Std.Dev.
     cluster (Intercept) 1.081    1.04   
    Number of obs: 6500, groups:  cluster, 849

    Fixed effects:
                              Estimate Std. Error z value Pr(>|z|)   
    (Intercept)               -2.29534    0.14281 -16.072   <2e-16 ***
    child_age                 -0.15959    0.07996  -1.996   0.0459 * 
    mother_educationprimary    0.12493    0.13413   0.931   0.3516   
    mother_educationsecondary  0.30202    0.15524   1.945   0.0517 . 
    mother_educationhigher     0.19440    0.33794   0.575   0.5651   
    —
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1

    Correlation of Fixed Effects:
                (Intr) chld_g mthr_dctnp mthr_dctns
    child_age   -0.284                            
    mthr_dctnpr -0.827  0.023                     
    mthr_dctnsc -0.732  0.009  0.759              
    mthr_dctnhg -0.340  0.013  0.346      0.323   <\/p>\n\n\n\n

    Keep models parsimonious<\/em>: Models should be kept as parsimonious as possible, including only predictors that are necessary and justified by theory or prior evidence. Complex models require more data to estimate parameters reliably, and overfitting can occur when too many variables are included relative to the number of outcome events. A parsimonious approach improves interpretability, reduces variance in estimates, and increases the likelihood that the model will generalise to other datasets or populations8<\/a><\/sup>.<\/p>\n\n\n\n

    Step 5: Model diagnostics and fit<\/em> Good modelling requires checking fit and influence<\/em>: Good modelling practice requires careful assessment of model fit and the influence of individual observations4<\/a><\/sup>. Multiple diagnostic tools should be employed rather than relying on a single test, as different diagnostics provide complementary information8<\/a><\/sup>. Evaluating residuals, leverage, influence measures, and overall goodness-of-fit help identify potential model misspecification, outliers, or highly influential points, ensuring that the final model is both robust and reliable.<\/p>\n\n\n\n

    Discrimination (AUC)<\/em>: Discrimination assesses a model\u2019s ability to distinguish between individuals who experience the event and those who do not. The c-statistic, commonly referred to as the area under the receiver operating characteristic curve (AUC), provides a summary measure of this separation50<\/a><\/sup>. In policy-relevant applications, an AUC above 0.7 is generally considered acceptable, while an AUC above 0.8 indicates strong discrimination; however, the interpretation should always consider the context, the outcome, and the intended use of the model51<\/a><\/sup>.<\/p>\n\n\n\n

    Calibration<\/em>: Calibration evaluates whether predicted probabilities align with observed outcomes52<\/a><\/sup>. Good calibration indicates that the model\u2019s predicted risk accurately reflects the actual probability of the event. Common approaches include calibration plots and the Hosmer-Lemeshow test46<\/a><\/sup>; the Hosmer-Lemeshow test can be overly sensitive in large samples, making visual inspection of calibration curves often more informative4<\/a><\/sup>. Properly calibrated models are critical for decision-making and policy applications, ensuring that predicted risks can be interpreted with confidence53<\/a><\/sup>.<\/p>\n\n\n\n

    R, Hosmer-Lemeshow, ROC, Brier:<\/p>\n\n\n\n

    library(ResourceSelection)<\/p>\n\n\n\n

    Warning: package ‘ResourceSelection’ was built under R version 4.5.1<\/p>\n\n\n\n

    ResourceSelection 0.3-6      2023-06-27<\/p>\n\n\n\n

    hoslem.test(dhs_mw$full_immunisation, fitted(mod_linear), g = 10)<\/p>\n\n\n\n


        Hosmer and Lemeshow goodness of fit (GOF) test

    data:  dhs_mw$full_immunisation, fitted(mod_linear)
    X-squared = 5.644, df = 8, p-value = 0.687<\/p>\n\n\n\n

    library(pROC)<\/p>\n\n\n\n

    Type ‘citation(“pROC”)’ for a citation.<\/p>\n\n\n\n


    Attaching package: ‘pROC’<\/p>\n\n\n\n

    The following objects are masked from ‘package:stats’:

        cov, smooth, var<\/p>\n\n\n\n

    roc_obj <- roc(dhs_mw$full_immunisation, fitted(mod_linear))<\/p>\n\n\n\n

    Setting levels: control = 0, case = 1<\/p>\n\n\n\n

    Setting direction: controls < cases<\/p>\n\n\n\n

    auc(roc_obj)<\/p>\n\n\n\n

    Area under the curve: 0.5484<\/p>\n\n\n\n

    # Brier score
    mean((dhs_mw$full_immunisation – fitted(mod_linear))^2)<\/p>\n\n\n\n

    1 0.1159879<\/p>\n\n\n\n

    Residuals and influential points<\/em>: Assessing residuals and influential points is essential to ensure model robustness4<\/a><\/sup>. Evaluating residuals, leverage, and influence measures helps identify observations that may disproportionately affect coefficient estimates55<\/a><\/sup>. This is particularly important in small samples, where single extreme or influential observations can distort results. Detecting and addressing such points through careful investigation or sensitivity analyses improves model reliability and strengthens confidence in the findings.<\/p>\n\n\n\n

    infl <- influence.measures(mod_linear)
    #summary(infl) #Activate this to see the results<\/p>\n\n\n\n

    Model validation<\/em>: Model validation is a critical step in both predictive and inferential analyses57<\/a><\/sup>. For prediction-focused models, internal validation can be performed using cross-validation or hold-out samples to assess performance and generalisability58<\/a><\/sup>. For inference-oriented models, bootstrapping provides a way to evaluate optimism, adjust confidence intervals, and check the stability of parameter estimates. In practice, tools such as the caret59<\/a><\/sup> or cv.glmnet60<\/a><\/sup> packages in R facilitate cross-validation for penalised models61<\/a><\/sup>, while the rms62<\/a><\/sup> package supports bootstrapping for inferential models. It is essential to report validation diagnostics in the manuscript, include relevant plots in supplementary materials, and document any influential observations that were removed or retained, ensuring transparency and reproducibility.<\/p>\n\n\n\n

    Step 6: Interpretation and presentation<\/em><\/p>\n\n\n\n

    Interpreting model results in a clear and actionable way is essential63<\/a><\/sup>, especially for policy and decision-making. While odds ratios are standard outputs in logistic regression, they can be difficult to interpret when the outcome is common64<\/a><\/sup>. Whenever possible, presenting predicted probabilities, absolute risks, or risk differences provides a more intuitive understanding of the effect sizes and facilitates practical decision-making65<\/a><\/sup>. Translating statistical findings into measures that are directly meaningful to stakeholders enhances the usability and impact of the research.<\/p>\n\n\n\n

    Odds ratios and confidence intervals, R:<\/p>\n\n\n\n

    library(broom)
    tidy(mod_linear, exponentiate = TRUE, conf.int = TRUE)<\/p>\n\n\n\n

    # A tibble: 10 \u00d7 7
       term                 estimate std.error statistic  p.value conf.low conf.high
       <chr>                   <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
     1 (Intercept)             0.132    0.131    -15.5   6.11e-54    0.102     0.170
     2 child_age               0.841    0.0731    -2.38  1.75e- 2    0.728     0.970
     3 mother_educationpri\u2026    1.04     0.121      0.344 7.31e- 1    0.826     1.33
     4 mother_educationsec\u2026    1.12     0.145      0.809 4.18e- 1    0.848     1.50
     5 mother_educationhig\u2026    0.967    0.308     -0.108 9.14e- 1    0.515     1.73
     6 wealth_quintilepoor\u2026    1.15     0.115      1.20  2.30e- 1    0.917     1.44
     7 wealth_quintilemidd\u2026    1.19     0.117      1.46  1.45e- 1    0.943     1.49
     8 wealth_quintilerich\u2026    1.37     0.119      2.63  8.55e- 3    1.08      1.73
     9 wealth_quintilerich\u2026    1.34     0.141      2.07  3.85e- 2    1.01      1.76
    10 urban1                  1.11     0.114      0.892 3.72e- 1    0.884     1.38<\/p>\n\n\n\n

    Predicted probabilities and marginal effects, R:<\/p>\n\n\n\n

    newdata <- expand.grid(child_age = c(1, 2), mother_education = c(“none”,”secondary”),
                           wealth_quintile = c(“poorest”,”richest”), urban = c(0,1))

    newdata$urban <- factor(newdata$urban, levels = levels(dhs_mw$urban))

    newdata$pred <- predict(mod_linear, newdata, type = “response”)
    newdata<\/p>\n\n\n\n

       child_age mother_education wealth_quintile urban       pred
    1          1             none         poorest     0 0.10006454
    2          2             none         poorest     0 0.08547725
    3          1        secondary         poorest     0 0.11117341
    4          2        secondary         poorest     0 0.09513788
    5          1             none         richest     0 0.12958323
    6          2             none         richest     0 0.11122458
    7          1        secondary         richest     0 0.14344683
    8          2        secondary         richest     0 0.12340252
    9          1             none         poorest     1 0.10959612
    10         2             none         poorest     1 0.09376410
    11         1        secondary         poorest     1 0.12162015
    12         2        secondary         poorest     1 0.10425454
    13         1             none         richest     1 0.14148453
    14         2             none         richest     1 0.12167547
    15         1        secondary         richest     1 0.15639262
    16         2        secondary         richest     1 0.13482409<\/p>\n\n\n\n

    Common outcomes<\/em>: When the outcome of interest is common, odds ratios can overstate the magnitude of associations, making interpretation challenging. In such cases, reporting adjusted risk ratios can provide a more intuitive measure of effect67<\/a><\/sup>. A practical approach is to use a Poisson regression model with robust standard errors4<\/a><\/sup>, which yields approximate risk ratios that are easier for policymakers and other stakeholders to interpret and act upon. This method improves the clarity and usability of findings without compromising statistical validity.<\/p>\n\n\n\n

    R, Poisson with robust SE:<\/p>\n\n\n\n

    library(sandwich)
    library(lmtest)<\/p>\n\n\n\n

    Loading required package: zoo<\/p>\n\n\n\n


    Attaching package: ‘zoo’<\/p>\n\n\n\n

    The following objects are masked from ‘package:base’:

        as.Date, as.Date.numeric<\/p>\n\n\n\n

    mod_pois <- glm(full_immunisation ~ child_age + mother_education + wealth_quintile + urban,
                    family = poisson(link = “log”), data = dhs_mw)
    coeftest(mod_pois, vcov = vcovHC(mod_pois, type = “HC0”))<\/p>\n\n\n\n


    z test of coefficients:

                               Estimate Std. Error  z value  Pr(>|z|)   
    (Intercept)               -2.148666   0.117840 -18.2337 < 2.2e-16 ***
    child_age                 -0.149931   0.063185  -2.3729  0.017649 * 
    mother_educationprimary    0.036564   0.105423   0.3468  0.728720   
    mother_educationsecondary  0.101486   0.125990   0.8055  0.420525   
    mother_educationhigher    -0.026540   0.263740  -0.1006  0.919846   
    wealth_quintilepoorer      0.120942   0.100385   1.2048  0.228285   
    wealth_quintilemiddle      0.149702   0.101941   1.4685  0.141963   
    wealth_quintilericher      0.271568   0.102762   2.6427  0.008225 **
    wealth_quintilerichest     0.253477   0.121070   2.0936  0.036292 * 
    urban1                     0.086229   0.094482   0.9126  0.361428   
    —
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1<\/p>\n\n\n\n

    Tables and figures<\/em>: Effective presentation of results is critical for clarity and impact. Tables should clearly display adjusted estimates, sample sizes, number of events, and key model fit statistics68<\/a><\/sup>. Supplementary tables can provide a complete list of variables considered, along with the code used to generate the final model, supporting transparency and reproducibility. Visualisations, such as predicted probability plots, are particularly valuable for communicating results to policymakers, translating complex statistical outputs into intuitive insights that can inform decision-making69<\/a><\/sup>.<\/p>\n\n\n\n

    Language and claims<\/em>: When reporting results from observational studies, careful attention should be paid to the language used70<\/a><\/sup>. Causal claims should be avoided unless the study design and methods explicitly support them71<\/a><\/sup>. Instead of stating that \u201cX causes Y,\u201d it is more appropriate to describe findings as \u201cX is associated with Y.\u201d For prediction-focused analyses, terms such as \u201cpredicts\u201d or \u201cforecast\u201d are suitable, reflecting the model\u2019s purpose without implying causation. Using cautious and precise language enhances credibility and prevents overinterpretation of the results17<\/a><\/sup>.<\/p>\n\n\n\n

    Practical checklist to finish modelling<\/strong><\/p>\n\n\n\n

    Before you write results, tick these items:<\/p>\n\n\n\n