I work with the risk of delays (0 = no risk, 1 = risk). I want to run a binary logistic regression considering panel data over time. I have data from 4 years (4 time-points) and some covariates (mother’s age — continuous; child’s gender — categorical,etc).
.
When I have just two time points, I can say: “Controlling for mother’s education and child`s age, the risk at the second time-point decreased or increased by “. However, with so many time-points, I don’t know if this interpretation is correct. So, I want to have an interpretation such as: “Controlling for mother education, child’s age (etc), the risk increased in 2020, 2021, etc”
I read about The glarma Package for Observation Driven Time Series Regression of Counts and the Package markovchain, but I could not figure out how to use them.
Code below:
library(tidyverse)
library(sjPlot)
df = structure(list(year_completed_cat = structure(c(3L, 4L, 4L, 3L,
3L, 2L, 4L, 2L, 2L, 2L, 2L, 3L, 4L, 4L, 2L, 5L, 3L, 3L, 4L, 2L,
5L, 2L, 3L, 3L, 3L, 4L, 2L, 3L, 2L, 3L, 2L, 4L, 3L, 4L, 4L, 2L,
2L, 5L, 4L, 3L, 3L, 2L, 3L, 5L, 5L, 3L, 3L, 5L, 2L, 3L, 2L, 2L,
2L, 3L, 4L, 3L, 2L, 2L, 2L, 2L, 2L, 3L, 2L, 3L, 2L, 2L, 2L, 2L,
3L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 2L, 3L, 3L, 5L, 3L, 5L, 3L, 3L,
3L, 2L, 5L, 3L, 3L, 4L, 2L, 4L, 3L, 3L, 3L, 5L, 2L, 2L, 2L, 3L,
3L, 3L, 2L, 3L, 3L, 2L, 3L, 2L, 3L, 3L, 3L, 4L, 3L, 3L, 3L, 2L,
2L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 4L, 5L, 2L, 3L, 5L, 3L, 3L, 5L,
2L, 2L, 4L, 2L, 3L, 2L, 5L, 3L, 2L, 2L, 2L, 3L, 2L, 4L, 2L, 3L,
3L, 5L), levels = c("18", "19", "20", "21", "22", "23", "24"), class = "factor"),
asqse_risk = c(0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1,
0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0,
1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0,
0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1,
1, 1, 0), momage = c(28, 28, 37, 23, 26, 20, 32, 32, 22,
32, 22, 32, 34, 31, 45, 29, 26, 24, 29, 27, 21, 32, 30, 27,
24, 27, 31, 36, 30, 33, 39, 25, 30, 27, 26, 23, 32, 33, 36,
36, 30, 29, 31, 29, 23, 22, 26, 29, 39, 19, 34, 33, 22, 32,
37, 30, 37, 30, 34, 27, 36, 19, 41, 35, 26, 28, 26, 27, 37,
29, 25, 25, 39, 31, 37, 26, 29, 33, 24, 24, 35, 22, 22, 23,
36, 28, 21, 22, 32, 34, 28, 34, 34, 28, 26, 34, 39, 31, 26,
22, 30, 35, 29, 32, 35, 28, 20, 28, 20, 34, 20, 24, 28, 22,
33, 36, 19, 24, 38, 26, 29, 29, 31, 25, 30, 25, 25, 21, 33,
29, 27, 30, 31, 30, 29, 38, 40, 33, 29, 28, 41, 30, 26, 23,
22, 31, 40, 26, 27, 18), sex_male = c(1, 1, 1, 0, 1, 1, 1,
0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0,
1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1,
1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1,
0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1,
1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1,
0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1,
1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0,
0, 0, 0, 0, 1, 1, 1, 0, 1, 0)), row.names = c(NA, -150L), class = "data.frame", na.action = structure(c(`4` = 4L,
`14` = 14L, `51` = 51L, `61` = 61L, `64` = 64L, `70` = 70L, `86` = 86L,
`109` = 109L, `113` = 113L, `115` = 115L, `156` = 156L, `159` = 159L,
`160` = 160L, `168` = 168L, `175` = 175L, `205` = 205L, `206` = 206L,
`221` = 221L, `222` = 222L, `227` = 227L, `235` = 235L, `241` = 241L,
`246` = 246L, `255` = 255L, `258` = 258L, `267` = 267L, `272` = 272L,
`273` = 273L, `276` = 276L, `298` = 298L, `304` = 304L, `311` = 311L,
`312` = 312L, `316` = 316L, `341` = 341L, `345` = 345L, `352` = 352L,
`375` = 375L, `388` = 388L, `391` = 391L, `393` = 393L, `400` = 400L,
`401` = 401L, `410` = 410L, `422` = 422L, `426` = 426L, `446` = 446L,
`450` = 450L, `492` = 492L, `504` = 504L, `506` = 506L, `520` = 520L,
`532` = 532L, `534` = 534L, `540` = 540L, `544` = 544L, `557` = 557L,
`565` = 565L, `595` = 595L, `609` = 609L, `612` = 612L, `622` = 622L,
`631` = 631L, `646` = 646L, `669` = 669L, `676` = 676L, `677` = 677L,
`685` = 685L, `710` = 710L, `722` = 722L, `725` = 725L, `730` = 730L,
`738` = 738L, `749` = 749L, `750` = 750L, `759` = 759L, `763` = 763L,
`768` = 768L, `774` = 774L, `777` = 777L, `781` = 781L, `784` = 784L,
`787` = 787L, `788` = 788L, `790` = 790L, `795` = 795L, `802` = 802L,
`803` = 803L, `807` = 807L, `816` = 816L, `819` = 819L, `833` = 833L,
`845` = 845L, `846` = 846L, `848` = 848L, `856` = 856L, `868` = 868L,
`871` = 871L, `891` = 891L, `915` = 915L, `919` = 919L, `921` = 921L,
`922` = 922L, `933` = 933L, `945` = 945L, `948` = 948L, `959` = 959L,
`961` = 961L, `967` = 967L, `969` = 969L, `970` = 970L, `974` = 974L,
`1000` = 1000L, `1007` = 1007L, `1018` = 1018L, `1022` = 1022L,
`1029` = 1029L, `1038` = 1038L, `1041` = 1041L, `1042` = 1042L,
`1061` = 1061L, `1070` = 1070L, `1075` = 1075L, `1081` = 1081L,
`1091` = 1091L, `1110` = 1110L, `1113` = 1113L, `1115` = 1115L,
`1131` = 1131L, `1137` = 1137L, `1139` = 1139L, `1144` = 1144L,
`1155` = 1155L, `1156` = 1156L, `1160` = 1160L, `1173` = 1173L,
`1178` = 1178L, `1182` = 1182L, `1185` = 1185L, `1188` = 1188L,
`1213` = 1213L, `1215` = 1215L, `1247` = 1247L, `1261` = 1261L,
`1267` = 1267L, `1273` = 1273L, `1279` = 1279L, `1282` = 1282L,
`1286` = 1286L, `1288` = 1288L, `1295` = 1295L, `1318` = 1318L,
`1335` = 1335L, `1354` = 1354L, `1370` = 1370L, `1372` = 1372L,
`1380` = 1380L, `1386` = 1386L, `1410` = 1410L, `1427` = 1427L,
`1443` = 1443L, `1453` = 1453L, `1461` = 1461L, `1488` = 1488L,
`1490` = 1490L, `1526` = 1526L, `1528` = 1528L, `1542` = 1542L,
`1573` = 1573L, `1578` = 1578L, `1596` = 1596L, `1612` = 1612L,
`1613` = 1613L, `1614` = 1614L, `1615` = 1615L, `1616` = 1616L,
`1617` = 1617L, `1618` = 1618L, `1619` = 1619L, `1628` = 1628L,
`1667` = 1667L, `1670` = 1670L, `1680` = 1680L, `1692` = 1692L,
`1705` = 1705L, `1733` = 1733L, `1741` = 1741L, `1751` = 1751L,
`1756` = 1756L, `1757` = 1757L, `1774` = 1774L, `1775` = 1775L,
`1778` = 1778L, `1785` = 1785L, `1793` = 1793L, `1799` = 1799L,
`1801` = 1801L, `1803` = 1803L, `1806` = 1806L, `1814` = 1814L,
`1825` = 1825L, `1828` = 1828L, `1831` = 1831L, `1832` = 1832L,
`1841` = 1841L, `1853` = 1853L, `1856` = 1856L, `1861` = 1861L,
`1868` = 1868L, `1872` = 1872L, `1874` = 1874L, `1896` = 1896L,
`1909` = 1909L, `1917` = 1917L, `1938` = 1938L, `1954` = 1954L,
`1963` = 1963L, `1964` = 1964L, `1971` = 1971L, `1994` = 1994L,
`1995` = 1995L, `2010` = 2010L, `2014` = 2014L, `2030` = 2030L,
`2034` = 2034L, `2036` = 2036L, `2043` = 2043L, `2051` = 2051L,
`2054` = 2054L, `2055` = 2055L, `2056` = 2056L, `2064` = 2064L,
`2066` = 2066L, `2068` = 2068L, `2089` = 2089L, `2093` = 2093L,
`2104` = 2104L, `2115` = 2115L), class = "omit"))
mod_longitudinal = df %>%
glm(asqse_risk ~ year_completed_cat ,
family = binomial(link = "logit"),
data = .)
sjPlot::tab_model(mod_longitudinal)
mod_longitudinal_covariates = df %>%
glm(asqse_risk ~ year_completed_cat + momage + sex_male,
family = binomial(link = "logit"),
data = .)
tab_model(mod_longitudinal_covariates)
tab_model(mod_longitudinal,mod_longitudinal_covariates)