Introduction to Latent Transition Analysis with R
This is an R Notebook, which means that I have prepared it using R Markdown so that you can read along and interact with the code in the boxes. When you conduct this yourself you will likely want to use RStudio, in that case the code in the boxes is used in the same way, just without all the extra explanations I have written for teaching purposes.
While we will not spend much time talking about how R works, instead we will talk about what you need to do specifically to run latent transition analysis. If you are interested you can find out more here: https://cran.r-project.org/doc/manuals/R-intro.pdf
Before we start you can see how r code works below.
First let’s ask r to do something simple, multiply 3 by 7. You can execute the code by pressing the green arrow or ctrl+shift+enter. If you want you can run individual lines by clicking on the line and pressing ctrl+shift.
3*7## [1] 21So we have our answer, but if we want to store it for later we can save it to an object. Here I have proposed “output” as the object and then asked r to store it with “<-”
output<-3*7That is now saved to “output.” We can check this by calling it.
output## [1] 21We can also now use that saved object in another calculation, here I want to add ten to the output:
output+10## [1] 31#You can also add notes that are not included in the code with "#"
#Packages have extra functions we can use.
#Below I demonstrate how to load packages, if it is not already installed then we need to install it first before we load it.
#Install package
if(!require(tidyverse)) +
install.packages("tidyverse")
#Load package
library(tidyverse)Tidyverse is an incredibly useful collection of packages that allow us to work with data in a tidy format. To find out more about tidy data you can read here: https://r4ds.hadley.nz/
Below I use a feature of the tidyverse package, the pipe “%>%”. This is an incredibly useful tool for chaining lines together. It essentially asks r to apply the outcome of the brackets in the line before the pipe to the brackets in the line after the pipe.
#Here I have created a function that squares the outcome of the brackets, I have assigned this to "square"
square <- function(x) x^2
#I then use the pipe to take the value stored in "output" and add 10 to it, in the next line we then apply the square function we defined earlier and apply it to the outcome
(output+10) %>%
square()## [1] 961Hopefully that all makes sense! We do not need to fully understand how R code works to run LTA, just what we need to change to run our own analyses.
Why Conduct LTA and When?
Setting Up
The first thing we want to do is install our packages, as we have already installed tidyverse we don’t need to do it again, however this code will check that for us.
Our central package is poLCA for the first portion of analysis, it allows us to conduct latent class analysis. While not strictly required for LTA we will use mice to help with missing data as it allows us to conduct multiple imputations. We will use readxl to read in the excel files that the data is stored in initially, flextable allows us to store dataframes as tables in a very readable format, and officer then lets us export tables cleanly to a word document.
#Install Packages
if(!require(poLCA)) +
install.packages("poLCA")
if(!require(mice)) +
install.packages("mice")
if(!require(readxl)) +
install.packages("readxl")
if(!require(officer)) +
install.packages("officer")
if(!require(flextable)) +
install.packages("flextable")Then we load the packages.
#Load Packages
library(tidyverse)
library(poLCA)
library(mice)
library(readxl)
library(officer)
library(flextable)Importing and Cleaning the Data
Data Importing
dataset_T1 <- read_excel("dataset_T1.xlsx")
dataset_T2 <- read_excel("dataset_T2.xlsx")Data Cleaning
#Inspecting the first dataset
names(dataset_T1)## [1] "Y1" "Y2" "Y3" "Y4" "Y5"
## [6] "Y6" "Y7" "Y8" "Y9" "Y10"
## [11] "deprivation" "gender" "score"summary(dataset_T1)## Y1 Y2 Y3 Y4
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:1.000
## Median :2.000 Median :2.000 Median :1.000 Median :2.000
## Mean :1.578 Mean :1.535 Mean :1.425 Mean :1.562
## 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:2.000
## Max. :2.000 Max. :2.000 Max. :2.000 Max. :2.000
## NA's :250 NA's :250
## Y5 Y6 Y7 Y8
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:1.000
## Median :2.000 Median :1.000 Median :2.000 Median :2.000
## Mean :1.562 Mean :1.391 Mean :1.575 Mean :1.578
## 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:2.000
## Max. :2.000 Max. :2.000 Max. :2.000 Max. :2.000
##
## Y9 Y10 deprivation gender
## Min. :1.000 Min. :1.000 Min. :1 Length:5000
## 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:2 Class :character
## Median :1.000 Median :1.000 Median :3 Mode :character
## Mean :1.428 Mean :1.415 Mean :3
## 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:4
## Max. :2.000 Max. :2.000 Max. :5
##
## score
## Length:5000
## Class :character
## Mode :character
##
##
##
## Lets loop a table for each column of the dataset, we will start with column 2 because 1 is ID.
for(i in 2:ncol(dataset_T1)) { #loops from column 2 onwards
print(table(dataset_T1[, i], useNA = "ifany")) #prints summary of all columns with NAs
cat("\n")} #adds a space## Y2
## 1 2 <NA>
## 2211 2539 250
##
## Y3
## 1 2
## 2876 2124
##
## Y4
## 1 2
## 2188 2812
##
## Y5
## 1 2
## 2192 2808
##
## Y6
## 1 2
## 3047 1953
##
## Y7
## 1 2
## 2126 2874
##
## Y8
## 1 2
## 2108 2892
##
## Y9
## 1 2
## 2858 2142
##
## Y10
## 1 2
## 2924 2076
##
## deprivation
## 1 2 3 4 5
## 1000 1000 1000 1000 1000
##
## gender
## 1 2
## 2500 2500
##
## score
## 1 2
## 2500 2500What do you notice?
Have a look at Y1_T1 and Y2_T1.
Now lets run it for the second dataset.
for(i in 2:ncol(dataset_T2)) {
print(table(dataset_T2[, i], useNA = "ifany"))
cat("\n")} ## Y2_T2
## 1 2 <NA>
## 2301 2449 250
##
## Y3_T2
## 1 2
## 2804 2196
##
## Y4_T2
## 1 2
## 2288 2712
##
## Y5_T2
## 1 2
## 2268 2732
##
## Y6_T2
## 1 2
## 2960 2040
##
## Y7_T2
## 1 2
## 2242 2758
##
## Y8_T2
## 1 2
## 2133 2867
##
## Y9_T2
## 1 2
## 2795 2205
##
## Y10_T2
## 1 2
## 2881 2119This dataset has the variables for T2 handily labelled differently from the previous dataset. As such we are safe to merge them before editing. Merging them is handy as it means we only need to keep track of one dataset.
We will use left_join from the tidyverse package to join the datasets. Left joining just means we are keeping everything from the first dataset, in this case dataset_T1 and adding what we choose from the second dataset. In this case we want the whole of the second dataset so the code is simple.
There is one other important consideration, we have unique identifiers for participants that we want to remain intact, so we need to tell R to join by that identifier, matching rows like for like. Unfortunately in the import the ID column is not stored in a way we can read.
We can verify that here:
head(dataset_T1)## # A tibble: 6 × 13
## Y1 Y2 Y3 Y4 Y5 Y6 Y7 Y8 Y9 Y10 deprivation gender
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 1 1 2 1 1 2 2 1 2 2 1 1
## 2 1 2 2 2 2 2 1 1 2 1 1 1
## 3 2 2 1 2 2 1 2 2 1 1 1 1
## 4 1 1 2 1 1 2 2 2 2 2 1 2
## 5 2 2 1 1 1 1 2 2 1 2 1 1
## 6 2 2 1 2 2 1 2 2 1 1 1 2
## # ℹ 1 more variable: score <chr>You can see that the first column is unnamed (or it may not even show up). So the first thing we will do is name it in a way that our code can see.
#We add the ID label here using mutate from the tidyverse package and relocate it to the first column (for aesthetic reasons)
dataset_T1 <- dataset_T1 %>%
mutate(ID = seq_len(nrow(.))) %>%
relocate(ID)
dataset_T2 <- dataset_T2 %>%
mutate(ID = seq_len(nrow(.))) %>%
relocate(ID)
#Then we left_join dataset_T1 and dataset_T2, by the ID variable we created
dataset_merged <- left_join(
dataset_T1,
dataset_T2,
by = "ID"
)Let’s check we got everything.
head(dataset_merged)## # A tibble: 6 × 24
## ID Y1 Y2 Y3 Y4 Y5 Y6 Y7 Y8 Y9 Y10 deprivation
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 2 1 1 2 2 1 2 2 1
## 2 2 1 2 2 2 2 2 1 1 2 1 1
## 3 3 2 2 1 2 2 1 2 2 1 1 1
## 4 4 1 1 2 1 1 2 2 2 2 2 1
## 5 5 2 2 1 1 1 1 2 2 1 2 1
## 6 6 2 2 1 2 2 1 2 2 1 1 1
## # ℹ 12 more variables: gender <chr>, score <chr>, Y1_T2 <dbl>, Y2_T2 <dbl>,
## # Y3_T2 <dbl>, Y4_T2 <dbl>, Y5_T2 <dbl>, Y6_T2 <dbl>, Y7_T2 <dbl>,
## # Y8_T2 <dbl>, Y9_T2 <dbl>, Y10_T2 <dbl>Now we are in a good place to clean this up in a way that stops us having to constantly look at our codebook!
#Let's add the codes for our factors with recode from the tidyverse package.
dataset_merged$gender <- recode(dataset_merged$gender,
"1" = "Male",
"2" = "Female")
dataset_merged$score <- recode(dataset_merged$score,
"1" = "Advanced",
"2" = "Beginner")
dataset_merged$deprivation <- recode(dataset_merged$deprivation,
"1" = "Quintile 1",
"2" = "Quintile 2",
"3" = "Quintile 3",
"4" = "Quintile 4",
"5" = "Quintile 5")Then lets check our work with the same loop from earlier. Again we ignore the first column as it contains ID.
for(i in 2:ncol(dataset_merged)) {
print(table(dataset_merged[, i], useNA = "ifany"))
cat("\n")}## Y1
## 1 2 <NA>
## 2006 2744 250
##
## Y2
## 1 2 <NA>
## 2211 2539 250
##
## Y3
## 1 2
## 2876 2124
##
## Y4
## 1 2
## 2188 2812
##
## Y5
## 1 2
## 2192 2808
##
## Y6
## 1 2
## 3047 1953
##
## Y7
## 1 2
## 2126 2874
##
## Y8
## 1 2
## 2108 2892
##
## Y9
## 1 2
## 2858 2142
##
## Y10
## 1 2
## 2924 2076
##
## deprivation
## Quintile 1 Quintile 2 Quintile 3 Quintile 4 Quintile 5
## 1000 1000 1000 1000 1000
##
## gender
## Female Male
## 2500 2500
##
## score
## Advanced Beginner
## 2500 2500
##
## Y1_T2
## 1 2 <NA>
## 2089 2661 250
##
## Y2_T2
## 1 2 <NA>
## 2301 2449 250
##
## Y3_T2
## 1 2
## 2804 2196
##
## Y4_T2
## 1 2
## 2288 2712
##
## Y5_T2
## 1 2
## 2268 2732
##
## Y6_T2
## 1 2
## 2960 2040
##
## Y7_T2
## 1 2
## 2242 2758
##
## Y8_T2
## 1 2
## 2133 2867
##
## Y9_T2
## 1 2
## 2795 2205
##
## Y10_T2
## 1 2
## 2881 2119We should see our labels being used now.
We have one more important piece of cleaning to do. You will note from the codebook that some questions are framed negatively and some are framed positively. We will want 1 and 2 to refer to the same positive or negative direction. As such we will flip the positive responses so that 1 = 2 and 2 = 1, while leaving the negative responses as they are. This should mean that 1 will always refer to the negative and 2 to the positive.
Before we do that though it will be helpful to rename the variables so we don’t have to keep working out which Y refers to what.
#Rename from the tidyverse is perfect for this
dataset_merged <- dataset_merged %>%
rename(
Tries_Best_T1 = Y1,
Enjoys_Learning_T1 = Y2,
Hates_school_T1 = Y3,
Looks_forward_T1 = Y4,
Feels_confident_T1 = Y5,
Gets_distracted_T1 = Y6,
Participates_T1 = Y7,
Homework_T1 = Y8,
Arrives_late_T1 = Y9,
Skips_School_T1 = Y10,
Tries_Best_T2 = Y1_T2,
Enjoys_Learning_T2 = Y2_T2,
Hates_school_T2 = Y3_T2,
Looks_forward_T2 = Y4_T2,
Feels_confident_T2 = Y5_T2,
Gets_distracted_T2 = Y6_T2,
Participates_T2 = Y7_T2,
Homework_T2 = Y8_T2,
Arrives_late_T2 = Y9_T2,
Skips_School_T2 = Y10_T2)
#quick check
names(dataset_merged)## [1] "ID" "Tries_Best_T1" "Enjoys_Learning_T1"
## [4] "Hates_school_T1" "Looks_forward_T1" "Feels_confident_T1"
## [7] "Gets_distracted_T1" "Participates_T1" "Homework_T1"
## [10] "Arrives_late_T1" "Skips_School_T1" "deprivation"
## [13] "gender" "score" "Tries_Best_T2"
## [16] "Enjoys_Learning_T2" "Hates_school_T2" "Looks_forward_T2"
## [19] "Feels_confident_T2" "Gets_distracted_T2" "Participates_T2"
## [22] "Homework_T2" "Arrives_late_T2" "Skips_School_T2"Now let’s flip the positive variables as discussed.
#We can use recode to do this
dataset_merged$Tries_Best_T1 <- recode(dataset_merged$Tries_Best_T1,
"1" = 2,
"2" = 1)
dataset_merged$Enjoys_Learning_T1 <- recode(dataset_merged$Enjoys_Learning_T1,
"1" = 2,
"2" = 1)
dataset_merged$Looks_forward_T1 <- recode(dataset_merged$Looks_forward_T1,
"1" = 2,
"2" = 1)
dataset_merged$Feels_confident_T1 <- recode(dataset_merged$Feels_confident_T1,
"1" = 2,
"2" = 1)
dataset_merged$Participates_T1 <- recode(dataset_merged$Participates_T1,
"1" = 2,
"2" = 1)
dataset_merged$Homework_T1 <- recode(dataset_merged$Homework_T1,
"1" = 2,
"2" = 1)
#And for T2
dataset_merged$Tries_Best_T2 <- recode(dataset_merged$Tries_Best_T2,
"1" = 2,
"2" = 1)
dataset_merged$Enjoys_Learning_T2 <- recode(dataset_merged$Enjoys_Learning_T2,
"1" = 2,
"2" = 1)
dataset_merged$Looks_forward_T2 <- recode(dataset_merged$Looks_forward_T2,
"1" = 2,
"2" = 1)
dataset_merged$Feels_confident_T2 <- recode(dataset_merged$Feels_confident_T2,
"1" = 2,
"2" = 1)
dataset_merged$Participates_T2 <- recode(dataset_merged$Participates_T2,
"1" = 2,
"2" = 1)
dataset_merged$Homework_T2 <- recode(dataset_merged$Homework_T2,
"1" = 2,
"2" = 1)Step 1. Latent Class Analysis
To conduct LTA, we will use the procedure designed by Lund and Ritz (2025) for use with R. Their 3 step process allows us to do what you used to need expensive software to do.
The first thing you would likely want to do is porduce some descriptive statistics for the dataset regarding responses to questions and overall characteristics. This can be done in any number of ways and is covered ad nauseam online and in other courses. We won’t go over that here but if you want to know more about how to do this in R you can read more with this brilliant resource from Danielle Navarro, 2019.
I also highly recommend checking out this resource on automating your tables by Leo Westbury, it is a little bit of extra effort to set up automated tables but once you start you will be wondering why you ever spent so much time copying over things by hand! We will be producing some automated tables as we go.
Now that our data is prepared we can conduct step 1: latent class analysis. We will need to do this twice, once at each time point. First though lets make lists of the variables we will include in the analysis.
Indicators_T1 <- cbind(Tries_Best_T1, Enjoys_Learning_T1, Hates_school_T1,
Looks_forward_T1, Feels_confident_T1, Gets_distracted_T1,
Participates_T1, Homework_T1, Arrives_late_T1, Skips_School_T1) ~ 1
Indicators_T2 <- cbind(Tries_Best_T2, Enjoys_Learning_T2, Hates_school_T2,
Looks_forward_T2, Feels_confident_T2, Gets_distracted_T2,
Participates_T2, Homework_T2, Arrives_late_T2, Skips_School_T2) ~ 1
#We need the ~ 1 to tell poLCA there are no predictors to account forWe should also set a seed to allow our analyses based on random number generation to be reproducible.
set.seed(968744)Then we run the analyses. We will run several, trying to fit an increasing number of classes to find the best fitting model. In this case we will investigate 2-5 latent engagement classes.
#After specifying we are using poLCA we then include the variable list first (Indicators_T1), then the dataset those variables can be found in (dataset_merged), then the number of latent classes we are trying to fit (nclass = 2), and finally the number of repetitions to attempt (nrep = 10). Finally we save the outcome to allow us to call it back later (lca_t1.2)
lca_t1.2 <- poLCA(Indicators_T1, dataset_merged, nclass = 2, nrep = 10)## Model 1: llik = -26458.78 ... best llik = -26458.78
## Model 2: llik = -26458.78 ... best llik = -26458.78
## Model 3: llik = -26458.78 ... best llik = -26458.78
## Model 4: llik = -26458.78 ... best llik = -26458.78
## Model 5: llik = -26458.78 ... best llik = -26458.78
## Model 6: llik = -26458.78 ... best llik = -26458.78
## Model 7: llik = -26458.78 ... best llik = -26458.78
## Model 8: llik = -26458.78 ... best llik = -26458.78
## Model 9: llik = -26458.78 ... best llik = -26458.78
## Model 10: llik = -26458.78 ... best llik = -26458.78
## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $Tries_Best_T1
## Pr(1) Pr(2)
## class 1: 0.1785 0.8215
## class 2: 0.8476 0.1524
##
## $Enjoys_Learning_T1
## Pr(1) Pr(2)
## class 1: 0.2176 0.7824
## class 2: 0.7475 0.2525
##
## $Hates_school_T1
## Pr(1) Pr(2)
## class 1: 0.1546 0.8454
## class 2: 0.8582 0.1418
##
## $Looks_forward_T1
## Pr(1) Pr(2)
## class 1: 0.3667 0.6333
## class 2: 0.6953 0.3047
##
## $Feels_confident_T1
## Pr(1) Pr(2)
## class 1: 0.2588 0.7412
## class 2: 0.7745 0.2255
##
## $Gets_distracted_T1
## Pr(1) Pr(2)
## class 1: 0.2523 0.7477
## class 2: 0.8525 0.1475
##
## $Participates_T1
## Pr(1) Pr(2)
## class 1: 0.3117 0.6883
## class 2: 0.7524 0.2476
##
## $Homework_T1
## Pr(1) Pr(2)
## class 1: 0.1704 0.8296
## class 2: 0.8617 0.1383
##
## $Arrives_late_T1
## Pr(1) Pr(2)
## class 1: 0.4190 0.5810
## class 2: 0.6702 0.3298
##
## $Skips_School_T1
## Pr(1) Pr(2)
## class 1: 0.2724 0.7276
## class 2: 0.7969 0.2031
##
## Estimated class population shares
## 0.4044 0.5956
##
## Predicted class memberships (by modal posterior prob.)
## 0.4031 0.5969
##
## =========================================================
## Fit for 2 latent classes:
## =========================================================
## number of observations: 4518
## number of estimated parameters: 21
## residual degrees of freedom: 1002
## maximum log-likelihood: -26458.78
##
## AIC(2): 52959.56
## BIC(2): 53094.29
## G^2(2): 1412.076 (Likelihood ratio/deviance statistic)
## X^2(2): 1492.148 (Chi-square goodness of fit)
## We will investigate this output later. For now we can be satisfied that no warning message occured. Now let’s do the same but for three latent engagement classes.
lca_t1.3 <- poLCA(Indicators_T1, dataset_merged, nclass = 3, nrep = 10)## Model 1: llik = -26339.08 ... best llik = -26339.08
## Model 2: llik = -26339.08 ... best llik = -26339.08
## Model 3: llik = -26339.08 ... best llik = -26339.08
## Model 4: llik = -26339.08 ... best llik = -26339.08
## Model 5: llik = -26339.08 ... best llik = -26339.08
## Model 6: llik = -26339.08 ... best llik = -26339.08
## Model 7: llik = -26339.08 ... best llik = -26339.08
## Model 8: llik = -26339.08 ... best llik = -26339.08
## Model 9: llik = -26339.08 ... best llik = -26339.08
## Model 10: llik = -26339.08 ... best llik = -26339.08
## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $Tries_Best_T1
## Pr(1) Pr(2)
## class 1: 0.8767 0.1233
## class 2: 0.5224 0.4776
## class 3: 0.0961 0.9039
##
## $Enjoys_Learning_T1
## Pr(1) Pr(2)
## class 1: 0.7895 0.2105
## class 2: 0.3997 0.6003
## class 3: 0.1814 0.8186
##
## $Hates_school_T1
## Pr(1) Pr(2)
## class 1: 0.8939 0.1061
## class 2: 0.4715 0.5285
## class 3: 0.0899 0.9101
##
## $Looks_forward_T1
## Pr(1) Pr(2)
## class 1: 0.6992 0.3008
## class 2: 0.6043 0.3957
## class 3: 0.2968 0.7032
##
## $Feels_confident_T1
## Pr(1) Pr(2)
## class 1: 0.8051 0.1949
## class 2: 0.4789 0.5211
## class 3: 0.2120 0.7880
##
## $Gets_distracted_T1
## Pr(1) Pr(2)
## class 1: 0.8944 0.1056
## class 2: 0.4770 0.5230
## class 3: 0.2086 0.7914
##
## $Participates_T1
## Pr(1) Pr(2)
## class 1: 0.7940 0.2060
## class 2: 0.4299 0.5701
## class 3: 0.2929 0.7071
##
## $Homework_T1
## Pr(1) Pr(2)
## class 1: 0.9014 0.0986
## class 2: 0.4829 0.5171
## class 3: 0.0977 0.9023
##
## $Arrives_late_T1
## Pr(1) Pr(2)
## class 1: 0.6846 0.3154
## class 2: 0.5375 0.4625
## class 3: 0.3893 0.6107
##
## $Skips_School_T1
## Pr(1) Pr(2)
## class 1: 0.8089 0.1911
## class 2: 0.6068 0.3932
## class 3: 0.1817 0.8183
##
## Estimated class population shares
## 0.5062 0.2013 0.2925
##
## Predicted class memberships (by modal posterior prob.)
## 0.5208 0.178 0.3012
##
## =========================================================
## Fit for 3 latent classes:
## =========================================================
## number of observations: 4518
## number of estimated parameters: 32
## residual degrees of freedom: 991
## maximum log-likelihood: -26339.08
##
## AIC(3): 52742.15
## BIC(3): 52947.46
## G^2(3): 1172.671 (Likelihood ratio/deviance statistic)
## X^2(3): 1054.569 (Chi-square goodness of fit)
## Again no error message, so let’s keep going.
lca_t1.4 <- poLCA(Indicators_T1, dataset_merged, nclass = 4, nrep = 10) # 4 classes## Model 1: llik = -26326.7 ... best llik = -26326.7
## Model 2: llik = -26328.59 ... best llik = -26326.7
## Model 3: llik = -26325.51 ... best llik = -26325.51
## Model 4: llik = -26327.63 ... best llik = -26325.51
## Model 5: llik = -26326.7 ... best llik = -26325.51
## Model 6: llik = -26328.27 ... best llik = -26325.51
## Model 7: llik = -26328.29 ... best llik = -26325.51
## Model 8: llik = -26328.37 ... best llik = -26325.51
## Model 9: llik = -26327.66 ... best llik = -26325.51
## Model 10: llik = -26329.55 ... best llik = -26325.51
## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $Tries_Best_T1
## Pr(1) Pr(2)
## class 1: 0.8768 0.1232
## class 2: 0.0966 0.9034
## class 3: 0.5160 0.4840
## class 4: 0.0900 0.9100
##
## $Enjoys_Learning_T1
## Pr(1) Pr(2)
## class 1: 0.7891 0.2109
## class 2: 0.2037 0.7963
## class 3: 0.3985 0.6015
## class 4: 0.1502 0.8498
##
## $Hates_school_T1
## Pr(1) Pr(2)
## class 1: 0.8944 0.1056
## class 2: 0.0753 0.9247
## class 3: 0.4597 0.5403
## class 4: 0.1092 0.8908
##
## $Looks_forward_T1
## Pr(1) Pr(2)
## class 1: 0.6994 0.3006
## class 2: 0.3055 0.6945
## class 3: 0.6023 0.3977
## class 4: 0.2796 0.7204
##
## $Feels_confident_T1
## Pr(1) Pr(2)
## class 1: 0.8044 0.1956
## class 2: 0.0000 1.0000
## class 3: 0.4806 0.5194
## class 4: 0.4476 0.5524
##
## $Gets_distracted_T1
## Pr(1) Pr(2)
## class 1: 0.8936 0.1064
## class 2: 0.2381 0.7619
## class 3: 0.4783 0.5217
## class 4: 0.1649 0.8351
##
## $Participates_T1
## Pr(1) Pr(2)
## class 1: 0.7936 0.2064
## class 2: 0.2288 0.7712
## class 3: 0.4275 0.5725
## class 4: 0.3663 0.6337
##
## $Homework_T1
## Pr(1) Pr(2)
## class 1: 0.9017 0.0983
## class 2: 0.0706 0.9294
## class 3: 0.4732 0.5268
## class 4: 0.1292 0.8708
##
## $Arrives_late_T1
## Pr(1) Pr(2)
## class 1: 0.6838 0.3162
## class 2: 0.4737 0.5263
## class 3: 0.5445 0.4555
## class 4: 0.2785 0.7215
##
## $Skips_School_T1
## Pr(1) Pr(2)
## class 1: 0.8092 0.1908
## class 2: 0.1981 0.8019
## class 3: 0.6012 0.3988
## class 4: 0.1561 0.8439
##
## Estimated class population shares
## 0.507 0.1548 0.2045 0.1337
##
## Predicted class memberships (by modal posterior prob.)
## 0.5224 0.2032 0.1868 0.0876
##
## =========================================================
## Fit for 4 latent classes:
## =========================================================
## number of observations: 4518
## number of estimated parameters: 43
## residual degrees of freedom: 980
## maximum log-likelihood: -26325.51
##
## AIC(4): 52737.02
## BIC(4): 53012.9
## G^2(4): 1145.539 (Likelihood ratio/deviance statistic)
## X^2(4): 1029.651 (Chi-square goodness of fit)
##
## ALERT: iterations finished, MAXIMUM LIKELIHOOD NOT FOUND
## We can see at the bottom that maximum likelihood was not found. To solve this we need to increase the maximum amount of iterations we allow the code to run. This will increase how long it takes but increases our chances of avoiding this error.
#We use maxiter to increase the maximum iterations
lca_t1.4 <- poLCA(Indicators_T1, dataset_merged, nclass = 4, nrep = 10, maxiter = 5000) # 4 classes## Model 1: llik = -26325.51 ... best llik = -26325.51
## Model 2: llik = -26326.7 ... best llik = -26325.51
## Model 3: llik = -26326.7 ... best llik = -26325.51
## Model 4: llik = -26328.15 ... best llik = -26325.51
## Model 5: llik = -26329.48 ... best llik = -26325.51
## Model 6: llik = -26326.7 ... best llik = -26325.51
## Model 7: llik = -26327.61 ... best llik = -26325.51
## Model 8: llik = -26326.7 ... best llik = -26325.51
## Model 9: llik = -26328.15 ... best llik = -26325.51
## Model 10: llik = -26326.7 ... best llik = -26325.51
## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $Tries_Best_T1
## Pr(1) Pr(2)
## class 1: 0.5160 0.4840
## class 2: 0.0967 0.9033
## class 3: 0.8768 0.1232
## class 4: 0.0900 0.9100
##
## $Enjoys_Learning_T1
## Pr(1) Pr(2)
## class 1: 0.3984 0.6016
## class 2: 0.2039 0.7961
## class 3: 0.7892 0.2108
## class 4: 0.1503 0.8497
##
## $Hates_school_T1
## Pr(1) Pr(2)
## class 1: 0.4597 0.5403
## class 2: 0.0753 0.9247
## class 3: 0.8944 0.1056
## class 4: 0.1090 0.8910
##
## $Looks_forward_T1
## Pr(1) Pr(2)
## class 1: 0.6023 0.3977
## class 2: 0.3057 0.6943
## class 3: 0.6994 0.3006
## class 4: 0.2795 0.7205
##
## $Feels_confident_T1
## Pr(1) Pr(2)
## class 1: 0.4807 0.5193
## class 2: 0.0000 1.0000
## class 3: 0.8044 0.1956
## class 4: 0.4455 0.5545
##
## $Gets_distracted_T1
## Pr(1) Pr(2)
## class 1: 0.4783 0.5217
## class 2: 0.2383 0.7617
## class 3: 0.8936 0.1064
## class 4: 0.1651 0.8349
##
## $Participates_T1
## Pr(1) Pr(2)
## class 1: 0.4276 0.5724
## class 2: 0.2283 0.7717
## class 3: 0.7935 0.2065
## class 4: 0.3662 0.6338
##
## $Homework_T1
## Pr(1) Pr(2)
## class 1: 0.4732 0.5268
## class 2: 0.0703 0.9297
## class 3: 0.9017 0.0983
## class 4: 0.1291 0.8709
##
## $Arrives_late_T1
## Pr(1) Pr(2)
## class 1: 0.5445 0.4555
## class 2: 0.4744 0.5256
## class 3: 0.6838 0.3162
## class 4: 0.2786 0.7214
##
## $Skips_School_T1
## Pr(1) Pr(2)
## class 1: 0.6012 0.3988
## class 2: 0.1984 0.8016
## class 3: 0.8092 0.1908
## class 4: 0.1560 0.8440
##
## Estimated class population shares
## 0.2045 0.1542 0.507 0.1343
##
## Predicted class memberships (by modal posterior prob.)
## 0.1864 0.2027 0.5224 0.0885
##
## =========================================================
## Fit for 4 latent classes:
## =========================================================
## number of observations: 4518
## number of estimated parameters: 43
## residual degrees of freedom: 980
## maximum log-likelihood: -26325.51
##
## AIC(4): 52737.02
## BIC(4): 53012.9
## G^2(4): 1145.538 (Likelihood ratio/deviance statistic)
## X^2(4): 1029.635 (Chi-square goodness of fit)
## #We will skip this because it will take too long but this is how we would normally run the 5 class model
#lca_t1.5 <- poLCA(Indicators_T1, dataset_merged, nclass = 5, nrep = 10, maxiter = 10000) # 5 classesTry it yourself! #1
Now try it yourself for time point 2. I will provide the skeleton, just fill in the correct Indicator list, dataset and class number.
#Here are the rest for Time point 2 as well
lca_t2.2 <- poLCA(Indicators_T2, dataset_merged, nclass = 2, nrep = 10) # 2 classes## Model 1: llik = -26742.56 ... best llik = -26742.56
## Model 2: llik = -26742.56 ... best llik = -26742.56
## Model 3: llik = -26742.56 ... best llik = -26742.56
## Model 4: llik = -26742.56 ... best llik = -26742.56
## Model 5: llik = -26742.56 ... best llik = -26742.56
## Model 6: llik = -26742.56 ... best llik = -26742.56
## Model 7: llik = -26742.56 ... best llik = -26742.56
## Model 8: llik = -26742.56 ... best llik = -26742.56
## Model 9: llik = -26742.56 ... best llik = -26742.56
## Model 10: llik = -26742.56 ... best llik = -26742.56
## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $Tries_Best_T2
## Pr(1) Pr(2)
## class 1: 0.1603 0.8397
## class 2: 0.8634 0.1366
##
## $Enjoys_Learning_T2
## Pr(1) Pr(2)
## class 1: 0.2250 0.7750
## class 2: 0.7381 0.2619
##
## $Hates_school_T2
## Pr(1) Pr(2)
## class 1: 0.1710 0.8290
## class 2: 0.8622 0.1378
##
## $Looks_forward_T2
## Pr(1) Pr(2)
## class 1: 0.3632 0.6368
## class 2: 0.6830 0.3170
##
## $Feels_confident_T2
## Pr(1) Pr(2)
## class 1: 0.2731 0.7269
## class 2: 0.7558 0.2442
##
## $Gets_distracted_T2
## Pr(1) Pr(2)
## class 1: 0.2574 0.7426
## class 2: 0.8457 0.1543
##
## $Participates_T2
## Pr(1) Pr(2)
## class 1: 0.3092 0.6908
## class 2: 0.7483 0.2517
##
## $Homework_T2
## Pr(1) Pr(2)
## class 1: 0.1922 0.8078
## class 2: 0.8646 0.1354
##
## $Arrives_late_T2
## Pr(1) Pr(2)
## class 1: 0.4037 0.5963
## class 2: 0.6772 0.3228
##
## $Skips_School_T2
## Pr(1) Pr(2)
## class 1: 0.2946 0.7054
## class 2: 0.7914 0.2086
##
## Estimated class population shares
## 0.4359 0.5641
##
## Predicted class memberships (by modal posterior prob.)
## 0.432 0.568
##
## =========================================================
## Fit for 2 latent classes:
## =========================================================
## number of observations: 4518
## number of estimated parameters: 21
## residual degrees of freedom: 1002
## maximum log-likelihood: -26742.56
##
## AIC(2): 53527.12
## BIC(2): 53661.85
## G^2(2): 1386.895 (Likelihood ratio/deviance statistic)
## X^2(2): 1424.194 (Chi-square goodness of fit)
## lca_t2.3 <- poLCA(Indicators_T2, dataset_merged, nclass = 3, nrep = 10, maxiter = 3000) # 3 classes## Model 1: llik = -26603.04 ... best llik = -26603.04
## Model 2: llik = -26603.04 ... best llik = -26603.04
## Model 3: llik = -26603.04 ... best llik = -26603.04
## Model 4: llik = -26603.04 ... best llik = -26603.04
## Model 5: llik = -26603.04 ... best llik = -26603.04
## Model 6: llik = -26603.04 ... best llik = -26603.04
## Model 7: llik = -26603.04 ... best llik = -26603.04
## Model 8: llik = -26603.04 ... best llik = -26603.04
## Model 9: llik = -26603.04 ... best llik = -26603.04
## Model 10: llik = -26603.04 ... best llik = -26603.04
## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $Tries_Best_T2
## Pr(1) Pr(2)
## class 1: 0.0918 0.9082
## class 2: 0.4691 0.5309
## class 3: 0.9066 0.0934
##
## $Enjoys_Learning_T2
## Pr(1) Pr(2)
## class 1: 0.1910 0.8090
## class 2: 0.3981 0.6019
## class 3: 0.7831 0.2169
##
## $Hates_school_T2
## Pr(1) Pr(2)
## class 1: 0.0935 0.9065
## class 2: 0.5035 0.4965
## class 3: 0.8980 0.1020
##
## $Looks_forward_T2
## Pr(1) Pr(2)
## class 1: 0.2830 0.7170
## class 2: 0.6129 0.3871
## class 3: 0.6846 0.3154
##
## $Feels_confident_T2
## Pr(1) Pr(2)
## class 1: 0.2214 0.7786
## class 2: 0.4890 0.5110
## class 3: 0.7868 0.2132
##
## $Gets_distracted_T2
## Pr(1) Pr(2)
## class 1: 0.2075 0.7925
## class 2: 0.4869 0.5131
## class 3: 0.8902 0.1098
##
## $Participates_T2
## Pr(1) Pr(2)
## class 1: 0.2967 0.7033
## class 2: 0.4136 0.5864
## class 3: 0.7961 0.2039
##
## $Homework_T2
## Pr(1) Pr(2)
## class 1: 0.1272 0.8728
## class 2: 0.4889 0.5111
## class 3: 0.9050 0.0950
##
## $Arrives_late_T2
## Pr(1) Pr(2)
## class 1: 0.3979 0.6021
## class 2: 0.4758 0.5242
## class 3: 0.7025 0.2975
##
## $Skips_School_T2
## Pr(1) Pr(2)
## class 1: 0.1841 0.8159
## class 2: 0.6509 0.3491
## class 3: 0.7992 0.2008
##
## Estimated class population shares
## 0.3123 0.2177 0.47
##
## Predicted class memberships (by modal posterior prob.)
## 0.3236 0.1923 0.4841
##
## =========================================================
## Fit for 3 latent classes:
## =========================================================
## number of observations: 4518
## number of estimated parameters: 32
## residual degrees of freedom: 991
## maximum log-likelihood: -26603.04
##
## AIC(3): 53270.08
## BIC(3): 53475.38
## G^2(3): 1107.849 (Likelihood ratio/deviance statistic)
## X^2(3): 999.1877 (Chi-square goodness of fit)
## lca_t2.4 <- poLCA(Indicators_T2, dataset_merged, nclass = 4, nrep = 10, maxiter = 6000) # 4 classes## Model 1: llik = -26591.14 ... best llik = -26591.14
## Model 2: llik = -26591.14 ... best llik = -26591.14
## Model 3: llik = -26589.8 ... best llik = -26589.8
## Model 4: llik = -26589.8 ... best llik = -26589.8
## Model 5: llik = -26589.8 ... best llik = -26589.8
## Model 6: llik = -26591.14 ... best llik = -26589.8
## Model 7: llik = -26589.8 ... best llik = -26589.8
## Model 8: llik = -26589.8 ... best llik = -26589.8
## Model 9: llik = -26595.83 ... best llik = -26589.8
## Model 10: llik = -26589.8 ... best llik = -26589.8
## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $Tries_Best_T2
## Pr(1) Pr(2)
## class 1: 0.7637 0.2363
## class 2: 0.9118 0.0882
## class 3: 0.0915 0.9085
## class 4: 0.4724 0.5276
##
## $Enjoys_Learning_T2
## Pr(1) Pr(2)
## class 1: 0.7440 0.2560
## class 2: 0.7841 0.2159
## class 3: 0.1908 0.8092
## class 4: 0.3999 0.6001
##
## $Hates_school_T2
## Pr(1) Pr(2)
## class 1: 0.9527 0.0473
## class 2: 0.8962 0.1038
## class 3: 0.0939 0.9061
## class 4: 0.5006 0.4994
##
## $Looks_forward_T2
## Pr(1) Pr(2)
## class 1: 1.0000 0.0000
## class 2: 0.6721 0.3279
## class 3: 0.2840 0.7160
## class 4: 0.6058 0.3942
##
## $Feels_confident_T2
## Pr(1) Pr(2)
## class 1: 0.0000 1.0000
## class 2: 0.8238 0.1762
## class 3: 0.2213 0.7787
## class 4: 0.4908 0.5092
##
## $Gets_distracted_T2
## Pr(1) Pr(2)
## class 1: 1.0000 0.0000
## class 2: 0.8868 0.1132
## class 3: 0.2083 0.7917
## class 4: 0.4815 0.5185
##
## $Participates_T2
## Pr(1) Pr(2)
## class 1: 0.9681 0.0319
## class 2: 0.7908 0.2092
## class 3: 0.2980 0.7020
## class 4: 0.4053 0.5947
##
## $Homework_T2
## Pr(1) Pr(2)
## class 1: 0.7991 0.2009
## class 2: 0.9090 0.0910
## class 3: 0.1268 0.8732
## class 4: 0.4912 0.5088
##
## $Arrives_late_T2
## Pr(1) Pr(2)
## class 1: 0.5812 0.4188
## class 2: 0.7076 0.2924
## class 3: 0.3978 0.6022
## class 4: 0.4771 0.5229
##
## $Skips_School_T2
## Pr(1) Pr(2)
## class 1: 0.6206 0.3794
## class 2: 0.8066 0.1934
## class 3: 0.1832 0.8168
## class 4: 0.6545 0.3455
##
## Estimated class population shares
## 0.0216 0.4486 0.3123 0.2175
##
## Predicted class memberships (by modal posterior prob.)
## 0.0139 0.4695 0.324 0.1926
##
## =========================================================
## Fit for 4 latent classes:
## =========================================================
## number of observations: 4518
## number of estimated parameters: 43
## residual degrees of freedom: 980
## maximum log-likelihood: -26589.8
##
## AIC(4): 53265.61
## BIC(4): 53541.49
## G^2(4): 1081.381 (Likelihood ratio/deviance statistic)
## X^2(4): 966.3725 (Chi-square goodness of fit)
## #Again we will skip the 5 class model for time reasons.
#lca_t2.5 <- poLCA(Indicators_T2, dataset_merged, nclass = 5, nrep = 10, maxiter = 10000) # 5 classesTo compare the fit of the models we will use a helper function that can harvest the relevant information from the outputs we stored.
#We save the function to ICfct to use later
ICfct <- function(lcafit) {
aicVal <- lcafit$aic
nparVal <- lcafit$npar
NobsVal <- lcafit$Nobs
llikVal <- lcafit$llik
aiccVal <- aicVal + (2 * nparVal^2 + 2 * nparVal) / (NobsVal - nparVal - 1)
abicVal <- nparVal * log((NobsVal + 2)/24) - 2 * llikVal
c(nons = NobsVal,
npar = nparVal,
llik = llikVal,
aic = aicVal,
aicc = aiccVal,
bic = lcafit$bic,
abic = abicVal)
}First let’s look at timepoint 1.
#Use the function on each T1 output - we could have added the fifth class here if we ran it
lc_t1.IC <- t(as.data.frame(lapply(list(lca_t1.2, lca_t1.3, lca_t1.4), ICfct)))
#Add row names 2-4 sequentially - adjust to suit the number of classes
row.names(lc_t1.IC) <- 2:4
#Round to 2 decimal places and print
round(lc_t1.IC, 2)## nons npar llik aic aicc bic abic
## 2 4518 21 -26458.78 52959.56 52959.76 53094.29 53027.56
## 3 4518 32 -26339.08 52742.15 52742.62 52947.46 52845.77
## 4 4518 43 -26325.51 52737.02 52737.86 53012.90 52876.26#Save the output as a dataframe
lc_t1.IC_df <- as.data.frame(lc_t1.IC)We are mostly interested in the fit statistics AIC and BIC, we also have adjusted versions. What we give most value to is complicated as discussed in the previous presentation. When looking at large datasets like this BIC is often preferred as it penalises for complexity, while AIC can increasingly predict large numbers of classes with larger datasets.
In this case AIC shows improvement between 2 and 3 classes, then slight improvement for each increasing number thereafter, however the difference is negligible (<10). BIC shows improvement for the 3 class solution then a drop off for 4 and 5 classes. In this case it seems reasonable to accept the 3 class solution.
Try it yourself! #2
Below, fill in the appropriate objects that we saved our lca outputs to. Remember we did not run the 5 class solution for time reasons.
lc_t2.IC <- t(as.data.frame(lapply(list(lca_t2.2, lca_t2.3, lca_t2.4), ICfct)))
row.names(lc_t2.IC) <- 2:4
round(lc_t2.IC, 2)## nons npar llik aic aicc bic abic
## 2 4518 21 -26742.56 53527.12 53527.33 53661.85 53595.12
## 3 4518 32 -26603.04 53270.08 53270.55 53475.38 53373.70
## 4 4518 43 -26589.80 53265.61 53266.45 53541.49 53404.85lc_t2.IC_df <- as.data.frame(lc_t2.IC)A similar picture emerges here, leaving us with a 3 class solution at each time point.
Sometimes this can be easier to see with visualisation. We will produce what are known as elbow plots. Visually we are looking for the point where increasing efficiency tails off (the elbow).
#For T1
#For our class labels
x_vals <- as.numeric(rownames(lc_t1.IC_df))
#This makes our plot labels universally horizontal
par(las = 1)
#We use matplot to quickly plot the output from our helper function
matplot(as.numeric(rownames(lc_t1.IC_df)), #here we specify the dataframe
lc_t1.IC_df[, c("aic", "bic")], #we plot AIC and BIC
type = "b", pch = 16:17, col = c("#1b9e77", "#e7298a"), lty = 1, lwd = 2, #aesthetic options
xlab = "Classes", ylab = "", xaxt = "n", main = "T1 Fit Statistics") #label names
axis(1, at = seq(min(x_vals), max(x_vals), by = 1)) #add x-axis with appropriate intervals
legend("right", legend = c("AIC","BIC"), #we also need the legend to interpret the plot
col = c("#1b9e77", "#e7298a"), pch = 16:17, lty = 1, lwd = 2, bty = "n")
The elbow becomes fairly obvious in this view.
Now let’s look at T2.
x_vals <- as.numeric(rownames(lc_t2.IC_df))
matplot(as.numeric(rownames(lc_t2.IC_df)),
lc_t2.IC_df[, c("aic", "bic")],
type = "b", pch = 16:17, col = c("#d95f02", "#7570b3"), lty = 1, lwd = 2,
xlab = "Classes", ylab = "", xaxt = "n", main = "T2 Fit Statistics")
axis(1, at = seq(min(x_vals), max(x_vals), by = 1))
legend("right", legend = c("AIC","BIC"),
col = c("#d95f02", "#7570b3"), pch = 16:17, lty = 1, lwd = 2, bty = "n")
Again we can see that elbow shape, especially for BIC.
Fit Statistic Table
The next thing we need to do is produce an output table for the fit statistics.
#You will recall that we created the datframe of fit stats already, so no need to redo this. We saved them as lc_t1.IC_df and lc_t2.IC_df.
#We will need to clean it up a bit with the column headers we want in order from left to right
colnames(lc_t1.IC_df) <- c("N", "Parameters", "Log-likelihood", "AIC", "AICc", "BIC", "aBIC")
colnames(lc_t2.IC_df) <- c("N", "Parameters", "Log-likelihood", "AIC", "AICc", "BIC", "aBIC")Now let’s make some flextables, these are very readable and clean.
#We will save this to ft_lca_t1, and tell R we are using the values stored in lca_t1.IC_df
ft_lca_t1 <- flextable(lc_t1.IC_df) %>%
theme_apa() %>% #This is a nice clean style
align(align = "center", part = "header") %>% #We will centre align
set_caption(caption = "Model Fit Statistics for T1", style = "italic") %>% #This is the title
autofit() #We fit the rows and columns to the size of the content
ft_lca_t1 #Then we view the tableN | Parameters | Log-likelihood | AIC | AICc | BIC | aBIC |
|---|---|---|---|---|---|---|
4,518.00 | 21.00 | -26,458.78 | 52,959.56 | 52,959.76 | 53,094.29 | 53,027.56 |
4,518.00 | 32.00 | -26,339.08 | 52,742.15 | 52,742.62 | 52,947.46 | 52,845.77 |
4,518.00 | 43.00 | -26,325.51 | 52,737.02 | 52,737.86 | 53,012.90 | 52,876.26 |
And again for time point 2
ft_lca_t2 <- flextable(lc_t2.IC_df) %>%
theme_apa() %>%
align(align = "center", part = "header") %>%
set_caption(caption = "Model Fit Statistics for T2", style = "italic") %>%
autofit()
ft_lca_t2N | Parameters | Log-likelihood | AIC | AICc | BIC | aBIC |
|---|---|---|---|---|---|---|
4,518.00 | 21.00 | -26,742.56 | 53,527.12 | 53,527.33 | 53,661.85 | 53,595.12 |
4,518.00 | 32.00 | -26,603.04 | 53,270.08 | 53,270.55 | 53,475.38 | 53,373.70 |
4,518.00 | 43.00 | -26,589.80 | 53,265.61 | 53,266.45 | 53,541.49 | 53,404.85 |
Finally we will use the officer package to export these to a single word document.
#This opens the process, telling officer to expect content to put in the doc.
Model_Fit_Tables <- read_docx()
#Then we add the tables we want
Model_Fit_Tables <- Model_Fit_Tables %>%
body_add_flextable(ft_lca_t1) %>% #add the table for T1
body_add_par("", style = "Normal") #we also add an empty paragraph to create a line space after each table
#Then add the next table
Model_Fit_Tables <- Model_Fit_Tables %>%
body_add_flextable(ft_lca_t2) %>%
body_add_par("", style = "Normal")
#And finally close the process, exporting this to the file name we specify in our working directory
print(Model_Fit_Tables, target = "Model Fit Tables.docx")Visualisation
Now we can investigate what our latent classes look like. We will use GGplot2 which allows us to highly customise what our outputs look like. But first we need the data in a useable format.
#First let's gather up all the probabilities for each variable from our LCA outputs
ir1 <- t(cbind(lca_t1.3[["probs"]][["Tries_Best_T1"]][, 2],
lca_t1.3[["probs"]][["Enjoys_Learning_T1"]][, 2],
lca_t1.3[["probs"]][["Hates_school_T1"]][, 2],
lca_t1.3[["probs"]][["Looks_forward_T1"]][, 2],
lca_t1.3[["probs"]][["Feels_confident_T1"]][, 2],
lca_t1.3[["probs"]][["Gets_distracted_T1"]][, 2],
lca_t1.3[["probs"]][["Participates_T1"]][, 2],
lca_t1.3[["probs"]][["Homework_T1"]][, 2],
lca_t1.3[["probs"]][["Arrives_late_T1"]][, 2],
lca_t1.3[["probs"]][["Skips_School_T1"]][, 2]))
ir2 <- t(cbind(lca_t2.3[["probs"]][["Tries_Best_T2"]][, 2],
lca_t2.3[["probs"]][["Enjoys_Learning_T2"]][, 2],
lca_t2.3[["probs"]][["Hates_school_T2"]][, 2],
lca_t2.3[["probs"]][["Looks_forward_T2"]][, 2],
lca_t2.3[["probs"]][["Feels_confident_T2"]][, 2],
lca_t2.3[["probs"]][["Gets_distracted_T2"]][, 2],
lca_t2.3[["probs"]][["Participates_T2"]][, 2],
lca_t2.3[["probs"]][["Homework_T2"]][, 2],
lca_t2.3[["probs"]][["Arrives_late_T2"]][, 2],
lca_t2.3[["probs"]][["Skips_School_T2"]][, 2]))
#Then we can order them any way we like, changing the numbers will change how they appear from left to right
new_order <- c(1,3,8,2,10,5,6,4,7,9)
#Lets also save the correct labels to apply to the graphs (\n means new line and will make these fit better on the graph)
indicator_labels <- c("Tries\nBest", "Enjoys\nLearning", "Hates\nschool", "Looks\nforward",
"Feels\nconfident", "Gets\ndistracted", "Participates",
"Homework", "Arrives\nlate", "Skips\nSchool")[new_order]
ir1_reordered <- ir1[new_order, ]
ir2_reordered <- ir2[new_order, ]Now we’re going to pivot our data longer, this means we are allowing the data to be grouped
ir1_reordered_tidy <- as.data.frame(ir1_reordered) %>%
mutate(Indicators = 1:10) %>% #First we create a column called indicators, labelled from 1-10
pivot_longer(
cols = -Indicators, #Then we pivot everything longer under indicators, so for each indicator there will be multiple rows, one for each latent class
names_to = "Group", #The latent class names will appear under Group
values_to = "Probability" #The actual proability values will appear under probability
)
#We will also add the labels we saved earlier in indicator_labels and tell R to treat the groups as factors
ir1_reordered_tidy <- ir1_reordered_tidy %>%
mutate(Indicators = factor(
Indicators,
levels = 1:10,
labels = indicator_labels))
ir1_reordered_tidy$Group <- factor(ir1_reordered_tidy$Group)
#Finally lets have a look
ir1_reordered_tidy## # A tibble: 30 × 3
## Indicators Group Probability
## <fct> <fct> <dbl>
## 1 "Tries\nBest" "class 1: " 0.123
## 2 "Tries\nBest" "class 2: " 0.478
## 3 "Tries\nBest" "class 3: " 0.904
## 4 "Hates\nschool" "class 1: " 0.106
## 5 "Hates\nschool" "class 2: " 0.528
## 6 "Hates\nschool" "class 3: " 0.910
## 7 "Homework" "class 1: " 0.0986
## 8 "Homework" "class 2: " 0.517
## 9 "Homework" "class 3: " 0.902
## 10 "Enjoys\nLearning" "class 1: " 0.210
## # ℹ 20 more rowsNow again for T2
ir2_reordered_tidy <- as.data.frame(ir2_reordered) %>%
mutate(Indicators = 1:10) %>%
pivot_longer(
cols = -Indicators,
names_to = "Group",
values_to = "Probability"
)
ir2_reordered_tidy <- ir2_reordered_tidy %>%
mutate(Indicators = factor(
Indicators,
levels = 1:10,
labels = indicator_labels))
ir2_reordered_tidy$Group <- factor(ir2_reordered_tidy$Group)
ir2_reordered_tidy## # A tibble: 30 × 3
## Indicators Group Probability
## <fct> <fct> <dbl>
## 1 "Tries\nBest" "class 1: " 0.908
## 2 "Tries\nBest" "class 2: " 0.531
## 3 "Tries\nBest" "class 3: " 0.0934
## 4 "Hates\nschool" "class 1: " 0.906
## 5 "Hates\nschool" "class 2: " 0.496
## 6 "Hates\nschool" "class 3: " 0.102
## 7 "Homework" "class 1: " 0.873
## 8 "Homework" "class 2: " 0.511
## 9 "Homework" "class 3: " 0.0950
## 10 "Enjoys\nLearning" "class 1: " 0.809
## # ℹ 20 more rowsNow one more dataset to make for a combined graph
#Add the identifier T1 for timepoint 1
ir1_df <- ir1_reordered_tidy %>%
mutate(Timepoint = "T1")
#Also rename the classes so we can label them, I found it quickest to trial and error this and you may want to choose the names once you see what the visualisation looks like!
ir1_df <- ir1_df %>%
mutate(Group = recode_factor(Group,
"class 3: " = "High",
"class 2: " = "Middle",
"class 1: " = "Low"))
#Now for T2
ir2_df <- ir2_reordered_tidy %>%
mutate(Timepoint = "T2")
ir2_df <- ir2_df %>%
mutate(Group = recode_factor(Group,
"class 1: " = "High",
"class 2: " = "Middle",
"class 3: " = "Low"))
#Then put them together
combined_df <- bind_rows(ir1_df, ir2_df)We are finally ready for visualisation!
First lets plot T1 with ggplot2 which is in the tidyverse package.
#First we name the plot, tell R we are using ggplot, then the dataframe we are using
#y is out y axis plot, and we are grouping by the variable we called Group
probplot_T1 <- ggplot(data=ir1_df, aes(x=Indicators, y=Probability, group=Group, shape = Group)) +
geom_line(colour = "#D81B60", linewidth = 1)+ #We choose the line colour and thickness
geom_point(colour = "#D81B60", size = 3)+ #And do the same for the individual points
ylim(0,1)+ #This is the scale for the y axis, from 0 to 1
labs(title = "Latent Groups For T1", x = "Indicator", y = "Item Response Probability")+ #labels
theme_classic()+ #A clean theme
theme(axis.text = element_text(size = 10))+ #font size of axis labels
scale_shape_manual(values = c(15,16,17)) #These are the point shapes - each number corresponds to a different shape
probplot_T1 #Then lets have a look
Try it yourself! #3
Now you try doing the same thing but for T2. You will need to plug in the correct df that we saved earlier.
#Then we do the same from T2
probplot_T2 <- ggplot(data=ir2_df, aes(x=Indicators, y=Probability, group=Group, shape = Group)) +
geom_line(colour = "#1E88E5", linewidth = 1)+
geom_point(colour = "#1E88E5", size = 3)+
ylim(0,1)+
labs(title = "Latent Groups For T2", x = "Indicator", y = "Item Response Probability")+
theme_classic()+
theme(axis.text = element_text(size = 10))+
scale_shape_manual(values = c(15,16,17))
probplot_T2
probplot_combined <- ggplot(data=combined_df, aes(x=Indicators, y=Probability, group=interaction(Timepoint, Group), colour = Timepoint, linetype = Timepoint, shape = Group)) + #here we need to tell R there is an interaction between time point and group we want to plot. We also differentiate colour and line type by timepoint and point shape by group
geom_line(linewidth = 1)+
geom_point(size = 3)+
ylim(0,1)+
labs(title = "Latent Groups For Primary Six and Secondary One", x = "Transitions Experiences Indicator", y = "Item Response Probability")+
theme_classic()+
theme(axis.text = element_text(size = 10))+
scale_colour_manual(values = c("T1" = "#D81B60", "T2" = "#1E88E5"))+ #assign a colour to each time point
scale_shape_manual(values = setNames(
c(15,16,17),
levels(combined_df$Group)))
probplot_combined #lets have a look
Other Tables
For our final part of step 1 we can a descriptive tables of how many are in each latent class. Classes are assigned based on maximum likelihood as we only actually have data on the likelihood of an individual being a member of the class.
I just want to quickly check which class is which. Because we have clear low, middle and high classes we can use the mean probabilities to work out which is which.
high_probs_1 <- sapply(lca_t1.3$probs, function(item) item[, 2])
class_score_1 <- rowMeans(high_probs_1)
class_score_1## class 1: class 2: class 3:
## 0.1852322 0.4989138 0.7953531high_probs_2 <- sapply(lca_t2.3$probs, function(item) item[, 2])
class_score_2 <- rowMeans(high_probs_2)
class_score_2## class 1: class 2: class 3:
## 0.7905772 0.5011293 0.1847990We can see that for T1 they are ordered low, middle, high but for T2 it’s high, middle, low.
We will create some dataframes to make the table but we will also use these as our primary categorical variables for the next steps of analysis.
#First we create an empty dataset with the same rows as dataset_merged
lcaT1 <- rep(NA, nrow(dataset_merged))
#Then we save our list of variables to bring over
complete_cases <- complete.cases(dataset_merged[, c(
"Tries_Best_T1", "Enjoys_Learning_T1", "Hates_school_T1", "Looks_forward_T1", "Feels_confident_T1",
"Gets_distracted_T1", "Participates_T1", "Homework_T1", "Arrives_late_T1", "Skips_School_T1")])
#Then for every participant that had complete data we bring over the predicted class, this preserves our missing data structure.
lcaT1[complete_cases] <- lca_t1.3$predclass
#Finally we assign each class number to the corresponding number from left to right and make those levels so they show in order - in this case 1:3 is right
lcaT1 <- factor(lcaT1, levels = 1:3, labels = c("Low", "Middle", "High"))
#Then do it again for T2
lcaT2 <- rep(NA, nrow(dataset_merged))
complete_cases <- complete.cases(dataset_merged[, c(
"Tries_Best_T2", "Enjoys_Learning_T2", "Hates_school_T2", "Looks_forward_T2", "Feels_confident_T2",
"Gets_distracted_T2", "Participates_T2", "Homework_T2", "Arrives_late_T2", "Skips_School_T2")])
lcaT2[complete_cases] <- lca_t2.3$predclass
lcaT2 <- factor(lcaT2, levels = 3:1, labels = c("Low", "Middle", "High")) #Here the order is reversed
#lets check if that looks about right
table(lcaT1)## lcaT1
## Low Middle High
## 2353 804 1361table(lcaT2)## lcaT2
## Low Middle High
## 2187 869 1462Now lets make the tables.
#Save the table form we just looked at as a dataframe
tab_T1 <- as.data.frame(table(lcaT1))
#Save the names for the top of our table columns
colnames(tab_T1) <- c("Class", "N")
#Calculate percentages, rounded to 2 decimal places
tab_T1$"Percentage (%)" <- round((tab_T1$N / sum(tab_T1$N)) * 100, 2)
#Save as an apa style flextable
ft_T1 <- flextable(tab_T1) %>% theme_apa() %>% fontsize(size = 11) %>% autofit()
#And view it
ft_T1Class | N | Percentage (%) |
|---|---|---|
Low | 2,353 | 52.08 |
Middle | 804 | 17.80 |
High | 1,361 | 30.12 |
And for T2.
tab_T2 <- as.data.frame(table(lcaT2))
colnames(tab_T2) <- c("Class", "N")
tab_T2$"Percentage (%)" <- round((tab_T2$N / sum(tab_T2$N)) * 100, 2)
ft_T2 <- flextable(tab_T2) %>% theme_apa() %>% fontsize(size = 11) %>% autofit()
ft_T2Class | N | Percentage (%) |
|---|---|---|
Low | 2,187 | 48.41 |
Middle | 869 | 19.23 |
High | 1,462 | 32.36 |
Then like before, save them together as a word document.
doc <- read_docx() %>%
body_add_par("Latent Class Membership by Time Point for T1", style = "heading 3") %>%
body_add_flextable(ft_T1) %>%
body_add_par("Latent Class Membership by Time Point for T2", style = "heading 3") %>%
body_add_flextable(ft_T2)
print(doc, target = "LTA Class Tables N.docx")Step 2. Transition Probabilities
Multiple Imputations
Before we begin step 2 we will approach some missing data. Because I simulated the dataset I know there is some missingness for two variables, however you may want to look into some resources for investigating missingness. E.g. Here. It is also worth noting that this is an adaptable design and if you have survey weights then you can apply them to the analysis. Be aware that you may need to look into alternative packages for conducting logistic regression with particularly complex survey weights.
We will use the pooled results of multiple imputed datasets to conduct our analyses. To do this we will use the mice package.
#First gather all the variables
ltadata.sim <- cbind(dataset_merged[, c("deprivation", "gender", "score",
"Tries_Best_T1", "Enjoys_Learning_T1", "Hates_school_T1", "Looks_forward_T1",
"Feels_confident_T1", "Gets_distracted_T1", "Participates_T1", "Homework_T1",
"Arrives_late_T1", "Skips_School_T1",
"Tries_Best_T2", "Enjoys_Learning_T2", "Hates_school_T2", "Looks_forward_T2",
"Feels_confident_T2", "Gets_distracted_T2", "Participates_T2", "Homework_T2",
"Arrives_late_T2", "Skips_School_T2")],
lcaT1, lcaT2) #We need to add these to the other variables in the dataframe#Just a quick check
str(ltadata.sim)## 'data.frame': 5000 obs. of 25 variables:
## $ deprivation : chr "Quintile 1" "Quintile 1" "Quintile 1" "Quintile 1" ...
## $ gender : chr "Male" "Male" "Male" "Female" ...
## $ score : chr "Advanced" "Beginner" "Advanced" "Beginner" ...
## $ Tries_Best_T1 : num 2 2 1 2 1 1 2 2 2 1 ...
## $ Enjoys_Learning_T1: num 2 1 1 2 1 1 2 2 2 1 ...
## $ Hates_school_T1 : num 2 2 1 2 1 1 2 2 2 1 ...
## $ Looks_forward_T1 : num 2 1 1 2 2 1 1 1 2 2 ...
## $ Feels_confident_T1: num 2 1 1 2 2 1 2 2 2 1 ...
## $ Gets_distracted_T1: num 2 2 1 2 1 1 2 2 1 1 ...
## $ Participates_T1 : num 1 2 1 1 1 1 2 1 1 1 ...
## $ Homework_T1 : num 2 2 1 1 1 1 2 2 2 1 ...
## $ Arrives_late_T1 : num 2 2 1 2 1 1 1 2 2 2 ...
## $ Skips_School_T1 : num 2 1 1 2 2 1 2 1 2 1 ...
## $ Tries_Best_T2 : num 2 2 1 2 1 1 2 2 2 2 ...
## $ Enjoys_Learning_T2: num 2 2 1 2 1 1 2 1 2 1 ...
## $ Hates_school_T2 : num 2 1 1 2 1 1 2 2 1 2 ...
## $ Looks_forward_T2 : num 1 2 1 2 1 1 2 2 1 1 ...
## $ Feels_confident_T2: num 1 2 1 2 1 1 1 1 2 2 ...
## $ Gets_distracted_T2: num 2 1 1 2 1 1 2 2 1 2 ...
## $ Participates_T2 : num 1 2 1 1 2 1 2 2 2 1 ...
## $ Homework_T2 : num 2 2 2 2 1 1 2 2 2 2 ...
## $ Arrives_late_T2 : num 1 1 1 1 1 2 2 2 2 1 ...
## $ Skips_School_T2 : num 2 1 1 2 1 2 2 1 2 2 ...
## $ lcaT1 : Factor w/ 3 levels "Low","Middle",..: 3 2 1 3 1 1 3 3 3 1 ...
## $ lcaT2 : Factor w/ 3 levels "Low","Middle",..: 3 2 1 3 1 1 3 3 3 3 ...Let’s make deprivation a factor with levels like gender and score.
ltadata.sim$deprivation <- factor(
ltadata.sim$deprivation,
levels = c("Quintile 1", "Quintile 2", "Quintile 3", "Quintile 4", "Quintile 5"),
labels = c("SIMD_1", "SIMD_2", "SIMD_3", "SIMD_4", "SIMD_5")
)
str(ltadata.sim)## 'data.frame': 5000 obs. of 25 variables:
## $ deprivation : Factor w/ 5 levels "SIMD_1","SIMD_2",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ gender : chr "Male" "Male" "Male" "Female" ...
## $ score : chr "Advanced" "Beginner" "Advanced" "Beginner" ...
## $ Tries_Best_T1 : num 2 2 1 2 1 1 2 2 2 1 ...
## $ Enjoys_Learning_T1: num 2 1 1 2 1 1 2 2 2 1 ...
## $ Hates_school_T1 : num 2 2 1 2 1 1 2 2 2 1 ...
## $ Looks_forward_T1 : num 2 1 1 2 2 1 1 1 2 2 ...
## $ Feels_confident_T1: num 2 1 1 2 2 1 2 2 2 1 ...
## $ Gets_distracted_T1: num 2 2 1 2 1 1 2 2 1 1 ...
## $ Participates_T1 : num 1 2 1 1 1 1 2 1 1 1 ...
## $ Homework_T1 : num 2 2 1 1 1 1 2 2 2 1 ...
## $ Arrives_late_T1 : num 2 2 1 2 1 1 1 2 2 2 ...
## $ Skips_School_T1 : num 2 1 1 2 2 1 2 1 2 1 ...
## $ Tries_Best_T2 : num 2 2 1 2 1 1 2 2 2 2 ...
## $ Enjoys_Learning_T2: num 2 2 1 2 1 1 2 1 2 1 ...
## $ Hates_school_T2 : num 2 1 1 2 1 1 2 2 1 2 ...
## $ Looks_forward_T2 : num 1 2 1 2 1 1 2 2 1 1 ...
## $ Feels_confident_T2: num 1 2 1 2 1 1 1 1 2 2 ...
## $ Gets_distracted_T2: num 2 1 1 2 1 1 2 2 1 2 ...
## $ Participates_T2 : num 1 2 1 1 2 1 2 2 2 1 ...
## $ Homework_T2 : num 2 2 2 2 1 1 2 2 2 2 ...
## $ Arrives_late_T2 : num 1 1 1 1 1 2 2 2 2 1 ...
## $ Skips_School_T2 : num 2 1 1 2 1 2 2 1 2 2 ...
## $ lcaT1 : Factor w/ 3 levels "Low","Middle",..: 3 2 1 3 1 1 3 3 3 1 ...
## $ lcaT2 : Factor w/ 3 levels "Low","Middle",..: 3 2 1 3 1 1 3 3 3 3 ...Now we can make the imputed datasets and save them to ltadata.imputed.
ltadata.imputed <- mice(ltadata.sim, m = 10) #m = 10 refers to the number of datasets we will produce. Higher degrees of missingness often require more datasets.##
## iter imp variable
## 1 1 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 1 2 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 1 3 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 1 4 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 1 5 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 1 6 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 1 7 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 1 8 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 1 9 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 1 10 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 2 1 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 2 2 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 2 3 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 2 4 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 2 5 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 2 6 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 2 7 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 2 8 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 2 9 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 2 10 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 3 1 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 3 2 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 3 3 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 3 4 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 3 5 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 3 6 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 3 7 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 3 8 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 3 9 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 3 10 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 4 1 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 4 2 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 4 3 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 4 4 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 4 5 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 4 6 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 4 7 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 4 8 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 4 9 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 4 10 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 5 1 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 5 2 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 5 3 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 5 4 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 5 5 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 5 6 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 5 7 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 5 8 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 5 9 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2
## 5 10 Tries_Best_T1 Enjoys_Learning_T1 Tries_Best_T2 Enjoys_Learning_T2 lcaT1 lcaT2## Warning: Number of logged events: 2Calculating Transition Probabilities
Now we will run a series of logistic regression analyses to calculate odds related to membership of a class at T2 in relation to class membership at T1. We will need to do this once per class at T2 as we will use each class as the outcome. If you have more than two time points then you can run this between each time point e.g. T1 -> T2; T2 -> T3, etc.
High class as outcome
We will run through the code for the T2 high class first, however the other analyses will look quite similar.
#First we will save this to ltadata.mice_high, then we need to use with because our dataset is a mids object from mice, as such we will be pooling the results of the analyses ran on each dataset.
#We then specify the name we have given the dataset and note we are conducting logisitc regression analysis by specifying glm and later family = binomial.
#Importantly we need to specify the outcome at T2 (lcaT2), in this case "High", and the predictors from T1 (lcaT1).
ltadata.mice_high <- with(data = ltadata.imputed, exp = glm(lcaT2 == "High" ~ lcaT1 - 1, family = binomial))
#Then we pool these analyses using Rubin's rules
logreg_high <- summary(pool(ltadata.mice_high))
#Let's have a look
logreg_high## term estimate std.error statistic df p.value
## 1 lcaT1Low -2.0801637 0.06320667 -32.910509 3664.146 2.447870e-208
## 2 lcaT1Middle -0.6103505 0.07188965 -8.490102 2008.731 3.951182e-17
## 3 lcaT1High 0.7045338 0.05530792 12.738388 4049.006 1.797624e-36We don’t need to interpret this yet, first we will transform this from a logarithmic scale to percentage odds as this is what is expected from LTA.
coef_high <- logreg_high[, 2] #Choosing the estimate column
llim_high <- coef_high - 1.96 * logreg_high[, 3] #Calculating the lower interval using standard error
ulim_high <- coef_high + 1.96 * logreg_high[, 3] #Calculating the upper interval
high_df <- data.frame(class = logreg_high[, 1],exp(cbind(coef_high, llim_high, ulim_high)) /
(1 + exp(cbind(coef_high, llim_high, ulim_high)))) #Pulling this together as a dataframe
high_df## class coef_high llim_high ulim_high
## 1 lcaT1Low 0.1110398 0.09938749 0.1238703
## 2 lcaT1Middle 0.3519793 0.32054808 0.3847472
## 3 lcaT1High 0.6691922 0.64477048 0.6927351We will look at interpreting this once we have done the same for the middle and lower groups.
Middle class as outcome
Try it yourself! #4
Here we are interested in “Middle” as the outcome, so try and plug it into the code.
ltadata.mice_mid <- with(data = ltadata.imputed, exp = glm(lcaT2 == "Middle" ~ lcaT1 - 1, family = binomial))
logreg_mid <- summary(pool(ltadata.mice_mid))
logreg_mid ## term estimate std.error statistic df p.value
## 1 lcaT1Low -1.8938478 0.05957982 -31.786735 2198.757 8.887925e-183
## 2 lcaT1Middle -0.4774431 0.06925465 -6.894022 4531.866 6.169905e-12
## 3 lcaT1High -1.4827430 0.06795294 -21.820144 2224.290 8.303889e-96coef_mid <- logreg_mid[, 2]
llim_mid <- coef_mid - 1.96 * logreg_mid[, 3]
ulim_mid <- coef_mid + 1.96 * logreg_mid[, 3]
mid_df <- data.frame(class = logreg_mid[, 1], exp(cbind(coef_mid, llim_mid, ulim_mid)) /
(1 + exp(cbind(coef_mid, llim_mid, ulim_mid))))
mid_df## class coef_mid llim_mid ulim_mid
## 1 lcaT1Low 0.1308064 0.1180919 0.1446651
## 2 lcaT1Middle 0.3828561 0.3513336 0.4153956
## 3 lcaT1High 0.1850135 0.1657668 0.2059431Low class as outcome
Try it yourself! #5
And do the same here but for “Low”.
ltadata.mice_low <- with(data = ltadata.imputed, exp = glm(lcaT2 == "Low" ~ lcaT1 - 1, family = binomial))
logreg_low <- summary(pool(ltadata.mice_low))
logreg_low## term estimate std.error statistic df p.value
## 1 lcaT1Low 1.142557 0.04683681 24.39442 2370.200 1.878640e-117
## 2 lcaT1Middle -1.019382 0.07721820 -13.20132 2783.729 1.219187e-38
## 3 lcaT1High -1.768053 0.07362816 -24.01328 4270.588 1.245106e-119coef_low <- logreg_low[, 2]
llim_low <- coef_low - 1.96 * logreg_low[, 3]
ulim_low <- coef_low + 1.96 * logreg_low[, 3]
low_df <- data.frame(class = logreg_low[, 1], exp(cbind(coef_low, llim_low, ulim_low)) /
(1 + exp(cbind(coef_low, llim_low, ulim_low))))
low_df## class coef_low llim_low ulim_low
## 1 lcaT1Low 0.7581487 0.7409201 0.7745802
## 2 lcaT1Middle 0.2651478 0.2367231 0.2956635
## 3 lcaT1High 0.1457846 0.1287154 0.1646894Make the Table
To interpret our step 2 findings, lets bring them together. We will be making another flextable
#This function can be used to make tables
make_table <- function(model_summary, T1_label) {
coef <- model_summary$estimate
se <- model_summary$std.error
llim <- coef - 1.96 * se
ulim <- coef + 1.96 * se
prob <- plogis(coef)
low <- plogis(llim)
high <- plogis(ulim)
data.frame(
T1 = T1_label,
T2 = model_summary$term,
prob = prob,
low = low,
high = high)}
#We use the function on each outcome
step2_df <- dplyr::bind_rows(
make_table(logreg_low, "Low"),
make_table(logreg_mid, "Middle"),
make_table(logreg_high, "High"))
#Quick recode so that the labels match
step2_df$T2 <- dplyr::recode(
step2_df$T2,
lcaT1Low = "Low",
lcaT1Middle = "Middle",
lcaT1High = "High")
#We need to multiply each estimate by 100 to get our percentages.
step2_df$label <- sprintf(
"%.1f%% (%.1f–%.1f)",
100 * step2_df$prob,
100 * step2_df$low,
100 * step2_df$high)
#Then we need to add out labels
table_wide <- step2_df %>%
dplyr::select(T1, T2, label) %>%
pivot_wider(names_from = T2, values_from = label)
table_wide <- table_wide[
match(c("Low", "Middle", "High"), table_wide$T1),
c("T1", "Low", "Middle", "High")]
#And make the flextable
step2_ft <- flextable(table_wide)
#Add some header labels to the table
step2_ft <- set_header_labels(
step2_ft,
T1 = "T1 Class",
Low = "T2 Low",
Middle = "T2 Middle",
High = "T2 High")
step2_ft <- autofit(step2_ft) #Fit the cells to the data
step2_ft #And have a lookT1 Class | T2 Low | T2 Middle | T2 High |
|---|---|---|---|
Low | 75.8% (74.1–77.5) | 26.5% (23.7–29.6) | 14.6% (12.9–16.5) |
Middle | 13.1% (11.8–14.5) | 38.3% (35.1–41.5) | 18.5% (16.6–20.6) |
High | 11.1% (9.9–12.4) | 35.2% (32.1–38.5) | 66.9% (64.5–69.3) |
You will note that the columns add to 100%, not the rows. As such you can see the relative odds of each outcome at T2, based on T1 class. We also have our 95% confidence intervals.
We can also save this to a word document.
doc <- read_docx() %>%
body_add_par("Transition Probabilities", style = "heading 3") %>%
body_add_flextable(step2_ft)
print(doc, target = "Transition Probabilities.docx")Step 3. Covariates
Our final step is adding covariates. This will allow us to look at transition probabilities from one class to another between time points, for a subgroup compared to a reference subgroup. For binary measures this is straightforward as there are only two groups to compare.
We run logistic regression models with each covariate for every transition we are interested in. Let’s start with gender, which in this case has been presented in binary.
## Outcome: Low to Middle
#First we come up with a way to save the data, it helps if we can keep track of this for making our tables.
#We need to use with because of our imputed dataset (ltadata.imputed)
#Then we use glm with family binomial for our logistic regression
#Importantly we input the later time point class first, so lcaT2 == "Middle"
#Then the timepoint one variable with the covariate, lcaT1*gender
#Lastly we ensure we look at a subset of the later time point to include the outcome of interest and also the same one as the earlier time point, this is because we are interested in odds of moving from low to middle versus staying low.
mice.step3.lomiG <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Middle" ~ lcaT1*gender, family = binomial,
subset = lcaT2 %in% c("Low", "Middle") & gender %in% c("Male", "Female")))We are interested only one line here, the rest can be disregarded for this method. We want the line that mentions the T1 class of interest (low in this case). The bottom three lines show the potential lines we could be interested in, we can see lcaT1Middle:genderMale and lcaT1High:genderMale, which are not of interest, in this case because low was the first level it is named simply genderMale (line 4). To allow us to interpret this, we must pool with Rubin’s rules and transform the data into the expected format. In this case we want odds ratios with 95% confidence intervals.
#Pool the results
logreg.step3.lomiG <- summary(pool(mice.step3.lomiG))
logreg.step3.lomiG## term estimate std.error statistic df
## 1 (Intercept) -1.58205357 0.07873784 -20.0926713 2820.696
## 2 lcaT1Middle 2.09278593 0.14623616 14.3110013 2535.052
## 3 lcaT1High 1.75862812 0.15032139 11.6991207 1656.148
## 4 genderMale -0.37898981 0.12148540 -3.1196326 1987.209
## 5 lcaT1Middle:genderMale 0.09551675 0.21025525 0.4542895 2459.888
## 6 lcaT1High:genderMale 0.50985531 0.22083703 2.3087401 2099.554
## p.value
## 1 4.941350e-84
## 2 9.849793e-45
## 3 1.956210e-30
## 4 1.836760e-03
## 5 6.496606e-01
## 6 2.105470e-02#Calculate the ORs from log odds
coefs.lomiG <- logreg.step3.lomiG
coefs.lomiG$OR <- exp(coefs.lomiG$estimate)
#Calculate the CIs
coefs.lomiG$OR_lower <- exp(coefs.lomiG$estimate - 1.96*coefs.lomiG$std.error)
coefs.lomiG$OR_upper <- exp(coefs.lomiG$estimate + 1.96*coefs.lomiG$std.error)
#Inspect the outcome
coefs.lomiG## term estimate std.error statistic df
## 1 (Intercept) -1.58205357 0.07873784 -20.0926713 2820.696
## 2 lcaT1Middle 2.09278593 0.14623616 14.3110013 2535.052
## 3 lcaT1High 1.75862812 0.15032139 11.6991207 1656.148
## 4 genderMale -0.37898981 0.12148540 -3.1196326 1987.209
## 5 lcaT1Middle:genderMale 0.09551675 0.21025525 0.4542895 2459.888
## 6 lcaT1High:genderMale 0.50985531 0.22083703 2.3087401 2099.554
## p.value OR OR_lower OR_upper
## 1 4.941350e-84 0.2055525 0.1761570 0.2398534
## 2 9.849793e-45 8.1074706 6.0870469 10.7985170
## 3 1.956210e-30 5.8044689 4.3232100 7.7932506
## 4 1.836760e-03 0.6845526 0.5395061 0.8685949
## 5 6.496606e-01 1.1002272 0.7286341 1.6613276
## 6 2.105470e-02 1.6650503 1.0800580 2.5668922We are interested only one line here, the rest can be disregarded for this method. We want the line that mentions the T1 class of interest (low in this case). The bottom three lines show the potential lines we could be interested in, we can see lcaT1Middle:genderMale and lcaT1High:genderMale, which are not of interest, in this case because low was the first level it is named simply genderMale (line 4). If we look at line 4 we can see the OR and CIs are 0.68 (0.54-0.87) meaning this is significant.
We can interpret this based on male being the reference. As the OR and CIs are all <1, this means that girls have significantly lower odds than boys of moving from the low to middle group, our findings suggest 32% lower odds.
We then do this for every transition of interest.
Try it yourself #6
Give it a go here for the low to high transition for gender. We need “High” from lcaT2, and we need to subset low and high from lcaT2.
## Outcome: low to high
mice.step3.lohiG <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "High" ~ lcaT1*gender, family = binomial,
subset = lcaT2 %in% c("Low", "High") & gender %in% c("Male", "Female")))
logreg.step3.lohiG <- summary(pool(mice.step3.lohiG))
logreg.step3.lohiG## term estimate std.error statistic df
## 1 (Intercept) -1.8611010 0.08933057 -20.8338654 2434.222
## 2 lcaT1Middle 2.2598710 0.15654992 14.4354653 1518.171
## 3 lcaT1High 3.2280799 0.13825395 23.3489158 1894.949
## 4 genderMale -0.1212994 0.12865513 -0.9428257 2353.689
## 5 lcaT1Middle:genderMale -0.1038859 0.21947547 -0.4733373 1755.154
## 6 lcaT1High:genderMale 0.4385044 0.20021517 2.1901657 2038.980
## p.value
## 1 7.695905e-89
## 2 2.324116e-44
## 3 3.482926e-106
## 4 3.458668e-01
## 5 6.360314e-01
## 6 2.862513e-02coefs.lohiG <- logreg.step3.lohiG
coefs.lohiG$OR <- exp(coefs.lohiG$estimate)
coefs.lohiG$OR_lower <- exp(coefs.lohiG$estimate - 1.96*coefs.lohiG$std.error)
coefs.lohiG$OR_upper <- exp(coefs.lohiG$estimate + 1.96*coefs.lohiG$std.error)
coefs.lohiG## term estimate std.error statistic df
## 1 (Intercept) -1.8611010 0.08933057 -20.8338654 2434.222
## 2 lcaT1Middle 2.2598710 0.15654992 14.4354653 1518.171
## 3 lcaT1High 3.2280799 0.13825395 23.3489158 1894.949
## 4 genderMale -0.1212994 0.12865513 -0.9428257 2353.689
## 5 lcaT1Middle:genderMale -0.1038859 0.21947547 -0.4733373 1755.154
## 6 lcaT1High:genderMale 0.4385044 0.20021517 2.1901657 2038.980
## p.value OR OR_lower OR_upper
## 1 7.695905e-89 0.1555013 0.1305252 0.1852567
## 2 2.324116e-44 9.5818530 7.0500390 13.0228935
## 3 3.482926e-106 25.2311645 19.2421304 33.0842608
## 4 3.458668e-01 0.8857688 0.6883462 1.1398135
## 5 6.360314e-01 0.9013281 0.5862213 1.3858117
## 6 2.862513e-02 1.5503867 1.0471610 2.2954437Gender
Here is the full code for every transition.
mice.step3.lomiG <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Middle" ~ lcaT1*gender, family = binomial,
subset = lcaT2 %in% c("Low", "Middle") & gender %in% c("Male", "Female")))
logreg.step3.lomiG <- summary(pool(mice.step3.lomiG))
logreg.step3.lomiG## term estimate std.error statistic df
## 1 (Intercept) -1.58205357 0.07873784 -20.0926713 2820.696
## 2 lcaT1Middle 2.09278593 0.14623616 14.3110013 2535.052
## 3 lcaT1High 1.75862812 0.15032139 11.6991207 1656.148
## 4 genderMale -0.37898981 0.12148540 -3.1196326 1987.209
## 5 lcaT1Middle:genderMale 0.09551675 0.21025525 0.4542895 2459.888
## 6 lcaT1High:genderMale 0.50985531 0.22083703 2.3087401 2099.554
## p.value
## 1 4.941350e-84
## 2 9.849793e-45
## 3 1.956210e-30
## 4 1.836760e-03
## 5 6.496606e-01
## 6 2.105470e-02coefs.lomiG <- logreg.step3.lomiG
coefs.lomiG$OR <- exp(coefs.lomiG$estimate)
coefs.lomiG$OR_lower <- exp(coefs.lomiG$estimate - 1.96*coefs.lomiG$std.error)
coefs.lomiG$OR_upper <- exp(coefs.lomiG$estimate + 1.96*coefs.lomiG$std.error)
coefs.lomiG## term estimate std.error statistic df
## 1 (Intercept) -1.58205357 0.07873784 -20.0926713 2820.696
## 2 lcaT1Middle 2.09278593 0.14623616 14.3110013 2535.052
## 3 lcaT1High 1.75862812 0.15032139 11.6991207 1656.148
## 4 genderMale -0.37898981 0.12148540 -3.1196326 1987.209
## 5 lcaT1Middle:genderMale 0.09551675 0.21025525 0.4542895 2459.888
## 6 lcaT1High:genderMale 0.50985531 0.22083703 2.3087401 2099.554
## p.value OR OR_lower OR_upper
## 1 4.941350e-84 0.2055525 0.1761570 0.2398534
## 2 9.849793e-45 8.1074706 6.0870469 10.7985170
## 3 1.956210e-30 5.8044689 4.3232100 7.7932506
## 4 1.836760e-03 0.6845526 0.5395061 0.8685949
## 5 6.496606e-01 1.1002272 0.7286341 1.6613276
## 6 2.105470e-02 1.6650503 1.0800580 2.5668922## Outcome: low to high
mice.step3.lohiG <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "High" ~ lcaT1*gender, family = binomial,
subset = lcaT2 %in% c("Low", "High") & gender %in% c("Male", "Female")))
logreg.step3.lohiG <- summary(pool(mice.step3.lohiG))
logreg.step3.lohiG## term estimate std.error statistic df
## 1 (Intercept) -1.8611010 0.08933057 -20.8338654 2434.222
## 2 lcaT1Middle 2.2598710 0.15654992 14.4354653 1518.171
## 3 lcaT1High 3.2280799 0.13825395 23.3489158 1894.949
## 4 genderMale -0.1212994 0.12865513 -0.9428257 2353.689
## 5 lcaT1Middle:genderMale -0.1038859 0.21947547 -0.4733373 1755.154
## 6 lcaT1High:genderMale 0.4385044 0.20021517 2.1901657 2038.980
## p.value
## 1 7.695905e-89
## 2 2.324116e-44
## 3 3.482926e-106
## 4 3.458668e-01
## 5 6.360314e-01
## 6 2.862513e-02coefs.lohiG <- logreg.step3.lohiG
coefs.lohiG$OR <- exp(coefs.lohiG$estimate)
coefs.lohiG$OR_lower <- exp(coefs.lohiG$estimate - 1.96*coefs.lohiG$std.error)
coefs.lohiG$OR_upper <- exp(coefs.lohiG$estimate + 1.96*coefs.lohiG$std.error)
coefs.lohiG## term estimate std.error statistic df
## 1 (Intercept) -1.8611010 0.08933057 -20.8338654 2434.222
## 2 lcaT1Middle 2.2598710 0.15654992 14.4354653 1518.171
## 3 lcaT1High 3.2280799 0.13825395 23.3489158 1894.949
## 4 genderMale -0.1212994 0.12865513 -0.9428257 2353.689
## 5 lcaT1Middle:genderMale -0.1038859 0.21947547 -0.4733373 1755.154
## 6 lcaT1High:genderMale 0.4385044 0.20021517 2.1901657 2038.980
## p.value OR OR_lower OR_upper
## 1 7.695905e-89 0.1555013 0.1305252 0.1852567
## 2 2.324116e-44 9.5818530 7.0500390 13.0228935
## 3 3.482926e-106 25.2311645 19.2421304 33.0842608
## 4 3.458668e-01 0.8857688 0.6883462 1.1398135
## 5 6.360314e-01 0.9013281 0.5862213 1.3858117
## 6 2.862513e-02 1.5503867 1.0471610 2.2954437## Outcome: middle to high
mice.step3.mihiG <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "High" ~ lcaT1*gender, family = binomial,
subset = lcaT2 %in% c("Middle", "High") & gender %in% c("Male", "Female")))
logreg.step3.mihiG <- summary(pool(mice.step3.mihiG))
logreg.step3.mihiG## term estimate std.error statistic df p.value
## 1 (Intercept) -0.2790474 0.1100198 -2.5363371 1854.222 1.128332e-02
## 2 lcaT1Middle 0.1670851 0.1563703 1.0685217 1597.759 2.854467e-01
## 3 lcaT1High 1.4694518 0.1467972 10.0100824 1888.315 5.099971e-23
## 4 genderMale 0.2576904 0.1651972 1.5598964 1515.505 1.189932e-01
## 5 lcaT1Middle:genderMale -0.1994027 0.2288755 -0.8712275 1801.388 3.837460e-01
## 6 lcaT1High:genderMale -0.0713509 0.2167586 -0.3291722 1309.721 7.420782e-01coefs.mihiG <- logreg.step3.mihiG
coefs.mihiG$OR <- exp(coefs.mihiG$estimate)
coefs.mihiG$OR_lower <- exp(coefs.mihiG$estimate - 1.96*coefs.mihiG$std.error)
coefs.mihiG$OR_upper <- exp(coefs.mihiG$estimate + 1.96*coefs.mihiG$std.error)
coefs.mihiG## term estimate std.error statistic df p.value
## 1 (Intercept) -0.2790474 0.1100198 -2.5363371 1854.222 1.128332e-02
## 2 lcaT1Middle 0.1670851 0.1563703 1.0685217 1597.759 2.854467e-01
## 3 lcaT1High 1.4694518 0.1467972 10.0100824 1888.315 5.099971e-23
## 4 genderMale 0.2576904 0.1651972 1.5598964 1515.505 1.189932e-01
## 5 lcaT1Middle:genderMale -0.1994027 0.2288755 -0.8712275 1801.388 3.837460e-01
## 6 lcaT1High:genderMale -0.0713509 0.2167586 -0.3291722 1309.721 7.420782e-01
## OR OR_lower OR_upper
## 1 0.7565040 0.6097622 0.938560
## 2 1.1818548 0.8698794 1.605718
## 3 4.3468515 3.2600070 5.796036
## 4 1.2939382 0.9360410 1.788678
## 5 0.8192200 0.5230916 1.282990
## 6 0.9311351 0.6088412 1.424037## Outcome: middle to low
mice.step3.miloG <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Low" ~ lcaT1*gender, family = binomial,
subset = lcaT2 %in% c("Middle", "Low") & gender %in% c("Male", "Female")))
logreg.step3.miloG <- summary(pool(mice.step3.miloG))
logreg.step3.miloG## term estimate std.error statistic df
## 1 (Intercept) 1.58205357 0.07873784 20.0926713 2820.696
## 2 lcaT1Middle -2.09278593 0.14623616 -14.3110013 2535.052
## 3 lcaT1High -1.75862812 0.15032139 -11.6991207 1656.148
## 4 genderMale 0.37898981 0.12148540 3.1196326 1987.209
## 5 lcaT1Middle:genderMale -0.09551675 0.21025525 -0.4542895 2459.888
## 6 lcaT1High:genderMale -0.50985531 0.22083703 -2.3087401 2099.554
## p.value
## 1 4.941350e-84
## 2 9.849793e-45
## 3 1.956210e-30
## 4 1.836760e-03
## 5 6.496606e-01
## 6 2.105470e-02coefs.miloG <- logreg.step3.miloG
coefs.miloG$OR <- exp(coefs.miloG$estimate)
coefs.miloG$OR_lower <- exp(coefs.miloG$estimate - 1.96*coefs.miloG$std.error)
coefs.miloG$OR_upper <- exp(coefs.miloG$estimate + 1.96*coefs.miloG$std.error)
coefs.miloG## term estimate std.error statistic df
## 1 (Intercept) 1.58205357 0.07873784 20.0926713 2820.696
## 2 lcaT1Middle -2.09278593 0.14623616 -14.3110013 2535.052
## 3 lcaT1High -1.75862812 0.15032139 -11.6991207 1656.148
## 4 genderMale 0.37898981 0.12148540 3.1196326 1987.209
## 5 lcaT1Middle:genderMale -0.09551675 0.21025525 -0.4542895 2459.888
## 6 lcaT1High:genderMale -0.50985531 0.22083703 -2.3087401 2099.554
## p.value OR OR_lower OR_upper
## 1 4.941350e-84 4.8649361 4.16921348 5.6767549
## 2 9.849793e-45 0.1233430 0.09260531 0.1642833
## 3 1.956210e-30 0.1722811 0.12831616 0.2313096
## 4 1.836760e-03 1.4608081 1.15128468 1.8535472
## 5 6.496606e-01 0.9089031 0.60192824 1.3724309
## 6 2.105470e-02 0.6005825 0.38957616 0.9258762## high to low
mice.step3.hiloG <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Low" ~ lcaT1*gender, family = binomial,
subset = lcaT2 %in% c("High", "Low") & gender %in% c("Male", "Female")))
logreg.step3.hiloG <- summary(pool(mice.step3.hiloG))
logreg.step3.hiloG## term estimate std.error statistic df
## 1 (Intercept) 1.8611010 0.08933057 20.8338654 2434.222
## 2 lcaT1Middle -2.2598710 0.15654992 -14.4354653 1518.171
## 3 lcaT1High -3.2280799 0.13825395 -23.3489158 1894.949
## 4 genderMale 0.1212994 0.12865513 0.9428257 2353.689
## 5 lcaT1Middle:genderMale 0.1038859 0.21947547 0.4733373 1755.154
## 6 lcaT1High:genderMale -0.4385044 0.20021517 -2.1901657 2038.980
## p.value
## 1 7.695905e-89
## 2 2.324116e-44
## 3 3.482926e-106
## 4 3.458668e-01
## 5 6.360314e-01
## 6 2.862513e-02coefs.hiloG <- logreg.step3.hiloG
coefs.hiloG$OR <- exp(coefs.hiloG$estimate)
coefs.hiloG$OR_lower <- exp(coefs.hiloG$estimate - 1.96*coefs.hiloG$std.error)
coefs.hiloG$OR_upper <- exp(coefs.hiloG$estimate + 1.96*coefs.hiloG$std.error)
coefs.hiloG## term estimate std.error statistic df
## 1 (Intercept) 1.8611010 0.08933057 20.8338654 2434.222
## 2 lcaT1Middle -2.2598710 0.15654992 -14.4354653 1518.171
## 3 lcaT1High -3.2280799 0.13825395 -23.3489158 1894.949
## 4 genderMale 0.1212994 0.12865513 0.9428257 2353.689
## 5 lcaT1Middle:genderMale 0.1038859 0.21947547 0.4733373 1755.154
## 6 lcaT1High:genderMale -0.4385044 0.20021517 -2.1901657 2038.980
## p.value OR OR_lower OR_upper
## 1 7.695905e-89 6.43081301 5.39791659 7.6613551
## 2 2.324116e-44 0.10436395 0.07678785 0.1418432
## 3 3.482926e-106 0.03963353 0.03022585 0.0519693
## 4 3.458668e-01 1.12896283 0.87733647 1.4527574
## 5 6.360314e-01 1.10947388 0.72159878 1.7058403
## 6 2.862513e-02 0.64500036 0.43564563 0.9549630## high to mid
mice.step3.himiG <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Middle" ~ lcaT1*gender, family = binomial,
subset = lcaT2 %in% c("High", "Middle") & gender %in% c("Male", "Female")))
logreg.step3.himiG <- summary(pool(mice.step3.himiG))
logreg.step3.himiG## term estimate std.error statistic df p.value
## 1 (Intercept) 0.2790474 0.1100198 2.5363371 1854.222 1.128332e-02
## 2 lcaT1Middle -0.1670851 0.1563703 -1.0685217 1597.759 2.854467e-01
## 3 lcaT1High -1.4694518 0.1467972 -10.0100824 1888.315 5.099971e-23
## 4 genderMale -0.2576904 0.1651972 -1.5598964 1515.505 1.189932e-01
## 5 lcaT1Middle:genderMale 0.1994027 0.2288755 0.8712275 1801.388 3.837460e-01
## 6 lcaT1High:genderMale 0.0713509 0.2167586 0.3291722 1309.721 7.420782e-01coefs.himiG <- logreg.step3.himiG
coefs.himiG$OR <- exp(coefs.himiG$estimate)
coefs.himiG$OR_lower <- exp(coefs.himiG$estimate - 1.96*coefs.himiG$std.error)
coefs.himiG$OR_upper <- exp(coefs.himiG$estimate + 1.96*coefs.himiG$std.error)
coefs.himiG## term estimate std.error statistic df p.value
## 1 (Intercept) 0.2790474 0.1100198 2.5363371 1854.222 1.128332e-02
## 2 lcaT1Middle -0.1670851 0.1563703 -1.0685217 1597.759 2.854467e-01
## 3 lcaT1High -1.4694518 0.1467972 -10.0100824 1888.315 5.099971e-23
## 4 genderMale -0.2576904 0.1651972 -1.5598964 1515.505 1.189932e-01
## 5 lcaT1Middle:genderMale 0.1994027 0.2288755 0.8712275 1801.388 3.837460e-01
## 6 lcaT1High:genderMale 0.0713509 0.2167586 0.3291722 1309.721 7.420782e-01
## OR OR_lower OR_upper
## 1 1.3218700 1.0654620 1.6399837
## 2 0.8461276 0.6227744 1.1495847
## 3 0.2300516 0.1725317 0.3067478
## 4 0.7728344 0.5590721 1.0683293
## 5 1.2206734 0.7794292 1.9117111
## 6 1.0739580 0.7022288 1.6424644Score
Then we can do the same for score, another binary covariate with only advanced and beginner to compare, so the code is identical with gender swapped for score and male/female swapped for advanced and beginner.
## Outcome: low to middle
mice.step3.lomiS <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Middle" ~ lcaT1*score, family = binomial,
subset = lcaT2 %in% c("Low", "Middle") & score %in% c("Advanced", "Beginner")))
logreg.step3.lomiS <- summary(pool(mice.step3.lomiS))
logreg.step3.lomiS## term estimate std.error statistic df
## 1 (Intercept) -1.77492895 0.0844381 -21.0204746 1717.040
## 2 lcaT1Middle 2.15400794 0.1481853 14.5359129 2401.199
## 3 lcaT1High 2.11216964 0.1561828 13.5236994 2275.916
## 4 scoreBeginner 0.03633059 0.1194582 0.3041281 2340.355
## 5 lcaT1Middle:scoreBeginner -0.05939340 0.2084400 -0.2849423 2699.549
## 6 lcaT1High:scoreBeginner -0.23022935 0.2185970 -1.0532137 2577.345
## p.value
## 1 1.763162e-87
## 2 6.048493e-46
## 3 3.855995e-40
## 4 7.610573e-01
## 5 7.757102e-01
## 6 2.923418e-01coefs.lomiS <- logreg.step3.lomiS
coefs.lomiS$OR <- exp(coefs.lomiS$estimate)
coefs.lomiS$OR_lower <- exp(coefs.lomiS$estimate - 1.96*coefs.lomiS$std.error)
coefs.lomiS$OR_upper <- exp(coefs.lomiS$estimate + 1.96*coefs.lomiS$std.error)
coefs.lomiS## term estimate std.error statistic df
## 1 (Intercept) -1.77492895 0.0844381 -21.0204746 1717.040
## 2 lcaT1Middle 2.15400794 0.1481853 14.5359129 2401.199
## 3 lcaT1High 2.11216964 0.1561828 13.5236994 2275.916
## 4 scoreBeginner 0.03633059 0.1194582 0.3041281 2340.355
## 5 lcaT1Middle:scoreBeginner -0.05939340 0.2084400 -0.2849423 2699.549
## 6 lcaT1High:scoreBeginner -0.23022935 0.2185970 -1.0532137 2577.345
## p.value OR OR_lower OR_upper
## 1 1.763162e-87 0.1694955 0.1436425 0.2000015
## 2 6.048493e-46 8.6193351 6.4466773 11.5242216
## 3 3.855995e-40 8.2661564 6.0863667 11.2266227
## 4 7.610573e-01 1.0369986 0.8205278 1.3105785
## 5 7.757102e-01 0.9423360 0.6262937 1.4178606
## 6 2.923418e-01 0.7943514 0.5175343 1.2192315## Outcome: low to high
mice.step3.lohiS <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "High" ~ lcaT1*score, family = binomial,
subset = lcaT2 %in% c("Low", "High") & score %in% c("Advanced", "Beginner")))
logreg.step3.lohiS <- summary(pool(mice.step3.lohiS))
logreg.step3.lohiS## term estimate std.error statistic df
## 1 (Intercept) -1.87666191 0.0873991 -21.4723247 2721.902
## 2 lcaT1Middle 2.18999547 0.1528399 14.3286934 2014.613
## 3 lcaT1High 3.38662849 0.1406563 24.0773399 3235.221
## 4 scoreBeginner -0.09404402 0.1279255 -0.7351468 3155.345
## 5 lcaT1Middle:scoreBeginner 0.03409506 0.2194745 0.1553486 1715.399
## 6 lcaT1High:scoreBeginner 0.12005671 0.1972999 0.6084987 3522.629
## p.value
## 1 1.300474e-94
## 2 2.042480e-44
## 3 5.765453e-118
## 4 4.623048e-01
## 5 8.765648e-01
## 6 5.428960e-01coefs.lohiS <- logreg.step3.lohiS
coefs.lohiS$OR <- exp(coefs.lohiS$estimate)
coefs.lohiS$OR_lower <- exp(coefs.lohiS$estimate - 1.96*coefs.lohiS$std.error)
coefs.lohiS$OR_upper <- exp(coefs.lohiS$estimate + 1.96*coefs.lohiS$std.error)
coefs.lohiS## term estimate std.error statistic df
## 1 (Intercept) -1.87666191 0.0873991 -21.4723247 2721.902
## 2 lcaT1Middle 2.18999547 0.1528399 14.3286934 2014.613
## 3 lcaT1High 3.38662849 0.1406563 24.0773399 3235.221
## 4 scoreBeginner -0.09404402 0.1279255 -0.7351468 3155.345
## 5 lcaT1Middle:scoreBeginner 0.03409506 0.2194745 0.1553486 1715.399
## 6 lcaT1High:scoreBeginner 0.12005671 0.1972999 0.6084987 3522.629
## p.value OR OR_lower OR_upper
## 1 1.300474e-94 0.1531003 0.1289973 0.181707
## 2 2.042480e-44 8.9351726 6.6222113 12.055990
## 3 5.765453e-118 29.5661017 22.4421802 38.951402
## 4 4.623048e-01 0.9102427 0.7083776 1.169633
## 5 8.765648e-01 1.0346830 0.6729562 1.590845
## 6 5.428960e-01 1.1275608 0.7659403 1.659912## Outcome: middle to high
mice.step3.mihiS <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "High" ~ lcaT1*score, family = binomial,
subset = lcaT2 %in% c("Middle", "High") & score %in% c("Advanced", "Beginner")))
logreg.step3.mihiS <- summary(pool(mice.step3.mihiS))
logreg.step3.mihiS## term estimate std.error statistic df
## 1 (Intercept) -0.10173296 0.1132375 -0.8984030 1510.577
## 2 lcaT1Middle 0.03598753 0.1600414 0.2248638 1503.794
## 3 lcaT1High 1.27445885 0.1480160 8.6102796 2128.887
## 4 scoreBeginner -0.13037462 0.1639573 -0.7951744 1620.012
## 5 lcaT1Middle:scoreBeginner 0.09348846 0.2309766 0.4047530 1217.697
## 6 lcaT1High:scoreBeginner 0.35028606 0.2120525 1.6518839 2268.184
## p.value
## 1 3.691139e-01
## 2 8.221158e-01
## 3 1.393093e-17
## 4 4.266286e-01
## 5 6.857302e-01
## 6 9.869664e-02coefs.mihiS <- logreg.step3.mihiS
coefs.mihiS$OR <- exp(coefs.mihiS$estimate)
coefs.mihiS$OR_lower <- exp(coefs.mihiS$estimate - 1.96*coefs.mihiS$std.error)
coefs.mihiS$OR_upper <- exp(coefs.mihiS$estimate + 1.96*coefs.mihiS$std.error)
coefs.mihiS## term estimate std.error statistic df
## 1 (Intercept) -0.10173296 0.1132375 -0.8984030 1510.577
## 2 lcaT1Middle 0.03598753 0.1600414 0.2248638 1503.794
## 3 lcaT1High 1.27445885 0.1480160 8.6102796 2128.887
## 4 scoreBeginner -0.13037462 0.1639573 -0.7951744 1620.012
## 5 lcaT1Middle:scoreBeginner 0.09348846 0.2309766 0.4047530 1217.697
## 6 lcaT1High:scoreBeginner 0.35028606 0.2120525 1.6518839 2268.184
## p.value OR OR_lower OR_upper
## 1 3.691139e-01 0.9032707 0.7234828 1.127737
## 2 8.221158e-01 1.0366429 0.7575288 1.418598
## 3 1.393093e-17 3.5767653 2.6760653 4.780620
## 4 4.266286e-01 0.8777665 0.6365255 1.210437
## 5 6.857302e-01 1.0979979 0.6982168 1.726684
## 6 9.869664e-02 1.4194735 0.9367519 2.150948## Outcome: middle to low
mice.step3.miloS <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Low" ~ lcaT1*score, family = binomial,
subset = lcaT2 %in% c("Middle", "Low") & score %in% c("Advanced", "Beginner")))
logreg.step3.miloS <- summary(pool(mice.step3.miloS))
logreg.step3.miloS## term estimate std.error statistic df
## 1 (Intercept) 1.77492895 0.0844381 21.0204746 1717.040
## 2 lcaT1Middle -2.15400794 0.1481853 -14.5359129 2401.199
## 3 lcaT1High -2.11216964 0.1561828 -13.5236994 2275.916
## 4 scoreBeginner -0.03633059 0.1194582 -0.3041281 2340.355
## 5 lcaT1Middle:scoreBeginner 0.05939340 0.2084400 0.2849423 2699.549
## 6 lcaT1High:scoreBeginner 0.23022935 0.2185970 1.0532137 2577.345
## p.value
## 1 1.763162e-87
## 2 6.048493e-46
## 3 3.855995e-40
## 4 7.610573e-01
## 5 7.757102e-01
## 6 2.923418e-01coefs.miloS <- logreg.step3.miloS
coefs.miloS$OR <- exp(coefs.miloS$estimate)
coefs.miloS$OR_lower <- exp(coefs.miloS$estimate - 1.96*coefs.miloS$std.error)
coefs.miloS$OR_upper <- exp(coefs.miloS$estimate + 1.96*coefs.miloS$std.error)
coefs.miloS## term estimate std.error statistic df
## 1 (Intercept) 1.77492895 0.0844381 21.0204746 1717.040
## 2 lcaT1Middle -2.15400794 0.1481853 -14.5359129 2401.199
## 3 lcaT1High -2.11216964 0.1561828 -13.5236994 2275.916
## 4 scoreBeginner -0.03633059 0.1194582 -0.3041281 2340.355
## 5 lcaT1Middle:scoreBeginner 0.05939340 0.2084400 0.2849423 2699.549
## 6 lcaT1High:scoreBeginner 0.23022935 0.2185970 1.0532137 2577.345
## p.value OR OR_lower OR_upper
## 1 1.763162e-87 5.8998619 4.99996180 6.9617274
## 2 6.048493e-46 0.1160182 0.08677376 0.1551187
## 3 3.855995e-40 0.1209752 0.08907398 0.1643016
## 4 7.610573e-01 0.9643214 0.76302181 1.2187277
## 5 7.757102e-01 1.0611926 0.70528797 1.5966950
## 6 2.923418e-01 1.2588887 0.82018876 1.9322391## high to low
mice.step3.hiloS <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Low" ~ lcaT1*score, family = binomial,
subset = lcaT2 %in% c("High", "Low") & score %in% c("Advanced", "Beginner")))
logreg.step3.hiloS <- summary(pool(mice.step3.hiloS))
logreg.step3.hiloS## term estimate std.error statistic df
## 1 (Intercept) 1.87666191 0.0873991 21.4723247 2721.902
## 2 lcaT1Middle -2.18999547 0.1528399 -14.3286934 2014.613
## 3 lcaT1High -3.38662849 0.1406563 -24.0773399 3235.221
## 4 scoreBeginner 0.09404402 0.1279255 0.7351468 3155.345
## 5 lcaT1Middle:scoreBeginner -0.03409506 0.2194745 -0.1553486 1715.399
## 6 lcaT1High:scoreBeginner -0.12005671 0.1972999 -0.6084987 3522.629
## p.value
## 1 1.300474e-94
## 2 2.042480e-44
## 3 5.765453e-118
## 4 4.623048e-01
## 5 8.765648e-01
## 6 5.428960e-01coefs.hiloS <- logreg.step3.hiloS
coefs.hiloS$OR <- exp(coefs.hiloS$estimate)
coefs.hiloS$OR_lower <- exp(coefs.hiloS$estimate - 1.96*coefs.hiloS$std.error)
coefs.hiloS$OR_upper <- exp(coefs.hiloS$estimate + 1.96*coefs.hiloS$std.error)
coefs.hiloS## term estimate std.error statistic df
## 1 (Intercept) 1.87666191 0.0873991 21.4723247 2721.902
## 2 lcaT1Middle -2.18999547 0.1528399 -14.3286934 2014.613
## 3 lcaT1High -3.38662849 0.1406563 -24.0773399 3235.221
## 4 scoreBeginner 0.09404402 0.1279255 0.7351468 3155.345
## 5 lcaT1Middle:scoreBeginner -0.03409506 0.2194745 -0.1553486 1715.399
## 6 lcaT1High:scoreBeginner -0.12005671 0.1972999 -0.6084987 3522.629
## p.value OR OR_lower OR_upper
## 1 1.300474e-94 6.53166512 5.50336468 7.75210289
## 2 2.042480e-44 0.11191726 0.08294632 0.15100696
## 3 5.765453e-118 0.03382252 0.02567302 0.04455895
## 4 4.623048e-01 1.09860811 0.85496908 1.41167651
## 5 8.765648e-01 0.96647962 0.62859685 1.48598082
## 6 5.428960e-01 0.88687014 0.60244165 1.30558477## high to mid
mice.step3.himiS <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Middle" ~ lcaT1*score, family = binomial,
subset = lcaT2 %in% c("High", "Middle") & score %in% c("Advanced", "Beginner")))
logreg.step3.himiS <- summary(pool(mice.step3.himiS))
logreg.step3.himiS## term estimate std.error statistic df
## 1 (Intercept) 0.10173296 0.1132375 0.8984030 1510.577
## 2 lcaT1Middle -0.03598753 0.1600414 -0.2248638 1503.794
## 3 lcaT1High -1.27445885 0.1480160 -8.6102796 2128.887
## 4 scoreBeginner 0.13037462 0.1639573 0.7951744 1620.012
## 5 lcaT1Middle:scoreBeginner -0.09348846 0.2309766 -0.4047530 1217.697
## 6 lcaT1High:scoreBeginner -0.35028606 0.2120525 -1.6518839 2268.184
## p.value
## 1 3.691139e-01
## 2 8.221158e-01
## 3 1.393093e-17
## 4 4.266286e-01
## 5 6.857302e-01
## 6 9.869664e-02coefs.himiS <- logreg.step3.himiS
coefs.himiS$OR <- exp(coefs.himiS$estimate)
coefs.himiS$OR_lower <- exp(coefs.himiS$estimate - 1.96*coefs.himiS$std.error)
coefs.himiS$OR_upper <- exp(coefs.himiS$estimate + 1.96*coefs.himiS$std.error)
coefs.himiS## term estimate std.error statistic df
## 1 (Intercept) 0.10173296 0.1132375 0.8984030 1510.577
## 2 lcaT1Middle -0.03598753 0.1600414 -0.2248638 1503.794
## 3 lcaT1High -1.27445885 0.1480160 -8.6102796 2128.887
## 4 scoreBeginner 0.13037462 0.1639573 0.7951744 1620.012
## 5 lcaT1Middle:scoreBeginner -0.09348846 0.2309766 -0.4047530 1217.697
## 6 lcaT1High:scoreBeginner -0.35028606 0.2120525 -1.6518839 2268.184
## p.value OR OR_lower OR_upper
## 1 3.691139e-01 1.1070878 0.8867319 1.382203
## 2 8.221158e-01 0.9646523 0.7049216 1.320082
## 3 1.393093e-17 0.2795822 0.2091779 0.373683
## 4 4.266286e-01 1.1392551 0.8261478 1.571029
## 5 6.857302e-01 0.9107485 0.5791449 1.432220
## 6 9.869664e-02 0.7044865 0.4649112 1.067518Deprivation Quintiles
In the case of deprivation there is a slight difference, it has 5 levels. As such we need a level to be the reference. This defaults to the first level which in this case is the least deprived quintile. When we were preparing the data earlier we could have reordered it if we wanted to compare to a different quintile. In this case comparing to the least deprived group seems reasonable so we don’t need to change anything. The difference will be that we are interested in multiple lines this time related to quintiles 2, 3, 4 and 5.
## Outcome: low to middle
mice.step3.lomiD <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Middle" ~ lcaT1*deprivation, family = binomial,
subset = lcaT2 %in% c("Low", "Middle")))
logreg.step3.lomiD <- summary(pool(mice.step3.lomiD))
logreg.step3.lomiD## term estimate std.error statistic df
## 1 (Intercept) -1.70284556 0.1311841 -12.9805777 2607.5331
## 2 lcaT1Middle 2.09204468 0.2399319 8.7193286 1375.5308
## 3 lcaT1High 1.83935211 0.2500233 7.3567215 1936.5442
## 4 deprivationSIMD_2 0.24964623 0.1778297 1.4038498 2817.0596
## 5 deprivationSIMD_3 0.06606173 0.1842109 0.3586200 2393.2834
## 6 deprivationSIMD_4 -0.54606173 0.2054523 -2.6578519 2840.5023
## 7 deprivationSIMD_5 -0.15083182 0.1891958 -0.7972261 2827.2780
## 8 lcaT1Middle:deprivationSIMD_2 -0.04428951 0.3330558 -0.1329793 1323.7855
## 9 lcaT1High:deprivationSIMD_2 0.03863502 0.3456650 0.1117701 2613.8558
## 10 lcaT1Middle:deprivationSIMD_3 -0.23163793 0.3470709 -0.6674080 919.4983
## 11 lcaT1High:deprivationSIMD_3 0.20886579 0.3487281 0.5989360 2350.9783
## 12 lcaT1Middle:deprivationSIMD_4 0.46373480 0.3410734 1.3596336 1507.2538
## 13 lcaT1High:deprivationSIMD_4 0.37944311 0.3579117 1.0601584 2531.8854
## 14 lcaT1Middle:deprivationSIMD_5 0.06889607 0.3358869 0.2051169 2102.0969
## 15 lcaT1High:deprivationSIMD_5 0.26420938 0.3541569 0.7460236 1127.5467
## p.value
## 1 2.212678e-37
## 2 7.918348e-18
## 3 2.769699e-13
## 4 1.604738e-01
## 5 7.199110e-01
## 6 7.908113e-03
## 7 4.253866e-01
## 8 8.942300e-01
## 9 9.110143e-01
## 10 5.046790e-01
## 11 5.492733e-01
## 12 1.741493e-01
## 13 2.891737e-01
## 14 8.375007e-01
## 15 4.558087e-01coefs.lomiD <- logreg.step3.lomiD
coefs.lomiD$OR <- exp(coefs.lomiD$estimate)
coefs.lomiD$OR_lower <- exp(coefs.lomiD$estimate - 1.96*coefs.lomiD$std.error)
coefs.lomiD$OR_upper <- exp(coefs.lomiD$estimate + 1.96*coefs.lomiD$std.error)
coefs.lomiD## term estimate std.error statistic df
## 1 (Intercept) -1.70284556 0.1311841 -12.9805777 2607.5331
## 2 lcaT1Middle 2.09204468 0.2399319 8.7193286 1375.5308
## 3 lcaT1High 1.83935211 0.2500233 7.3567215 1936.5442
## 4 deprivationSIMD_2 0.24964623 0.1778297 1.4038498 2817.0596
## 5 deprivationSIMD_3 0.06606173 0.1842109 0.3586200 2393.2834
## 6 deprivationSIMD_4 -0.54606173 0.2054523 -2.6578519 2840.5023
## 7 deprivationSIMD_5 -0.15083182 0.1891958 -0.7972261 2827.2780
## 8 lcaT1Middle:deprivationSIMD_2 -0.04428951 0.3330558 -0.1329793 1323.7855
## 9 lcaT1High:deprivationSIMD_2 0.03863502 0.3456650 0.1117701 2613.8558
## 10 lcaT1Middle:deprivationSIMD_3 -0.23163793 0.3470709 -0.6674080 919.4983
## 11 lcaT1High:deprivationSIMD_3 0.20886579 0.3487281 0.5989360 2350.9783
## 12 lcaT1Middle:deprivationSIMD_4 0.46373480 0.3410734 1.3596336 1507.2538
## 13 lcaT1High:deprivationSIMD_4 0.37944311 0.3579117 1.0601584 2531.8854
## 14 lcaT1Middle:deprivationSIMD_5 0.06889607 0.3358869 0.2051169 2102.0969
## 15 lcaT1High:deprivationSIMD_5 0.26420938 0.3541569 0.7460236 1127.5467
## p.value OR OR_lower OR_upper
## 1 2.212678e-37 0.1821644 0.1408632 0.2355753
## 2 7.918348e-18 8.1014631 5.0620839 12.9657480
## 3 2.769699e-13 6.2924601 3.8547508 10.2717546
## 4 1.604738e-01 1.2835712 0.9058332 1.8188283
## 5 7.199110e-01 1.0682927 0.7445378 1.5328292
## 6 7.908113e-03 0.5792265 0.3872255 0.8664289
## 7 4.253866e-01 0.8599923 0.5935371 1.2460666
## 8 8.942300e-01 0.9566770 0.4980381 1.8376724
## 9 9.110143e-01 1.0393911 0.5278894 2.0465152
## 10 5.046790e-01 0.7932333 0.4017614 1.5661511
## 11 5.492733e-01 1.2322796 0.6221082 2.4409145
## 12 1.741493e-01 1.5900013 0.8148355 3.1025942
## 13 2.891737e-01 1.4614705 0.7246519 2.9474785
## 14 8.375007e-01 1.0713249 0.5546365 2.0693497
## 15 4.558087e-01 1.3024009 0.6505493 2.6074087## Outcome: low to high
mice.step3.lohiD <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "High" ~ lcaT1*deprivation, family = binomial,
subset = lcaT2 %in% c("Low", "High")))
logreg.step3.lohiD <- summary(pool(mice.step3.lohiD))
logreg.step3.lohiD## term estimate std.error statistic df
## 1 (Intercept) -1.74625978 0.1333497 -13.09533868 3201.180
## 2 lcaT1Middle 2.05223339 0.2410168 8.51489685 2486.215
## 3 lcaT1High 3.31298684 0.2145886 15.43878319 3393.767
## 4 deprivationSIMD_2 0.01092666 0.1893618 0.05770256 3585.001
## 5 deprivationSIMD_3 0.07106083 0.1863199 0.38139143 3333.752
## 6 deprivationSIMD_4 -0.64834552 0.2149082 -3.01684933 3534.271
## 7 deprivationSIMD_5 -0.44423859 0.2078080 -2.13773583 2921.773
## 8 lcaT1Middle:deprivationSIMD_2 0.12866793 0.3408672 0.37747228 2340.421
## 9 lcaT1High:deprivationSIMD_2 -0.01633763 0.3082935 -0.05299377 3786.868
## 10 lcaT1Middle:deprivationSIMD_3 -0.12505026 0.3409031 -0.36682051 2661.730
## 11 lcaT1High:deprivationSIMD_3 0.04688613 0.3061281 0.15315854 3174.414
## 12 lcaT1Middle:deprivationSIMD_4 0.46230947 0.3470174 1.33223722 3316.506
## 13 lcaT1High:deprivationSIMD_4 0.43194974 0.3145081 1.37341372 3852.284
## 14 lcaT1Middle:deprivationSIMD_5 0.44409275 0.3478262 1.27676626 2524.584
## 15 lcaT1High:deprivationSIMD_5 0.35219871 0.3163132 1.11344922 2393.207
## p.value
## 1 3.301268e-38
## 2 2.840762e-17
## 3 5.065286e-52
## 4 9.539888e-01
## 5 7.029371e-01
## 6 2.572364e-03
## 7 3.262095e-02
## 8 7.058569e-01
## 9 9.577397e-01
## 10 7.137821e-01
## 11 8.782830e-01
## 12 1.828738e-01
## 13 1.697037e-01
## 14 2.018022e-01
## 15 2.656274e-01coefs.lohiD <- logreg.step3.lohiD
coefs.lohiD$OR <- exp(coefs.lohiD$estimate)
coefs.lohiD$OR_lower <- exp(coefs.lohiD$estimate - 1.96*coefs.lohiD$std.error)
coefs.lohiD$OR_upper <- exp(coefs.lohiD$estimate + 1.96*coefs.lohiD$std.error)
coefs.lohiD## term estimate std.error statistic df
## 1 (Intercept) -1.74625978 0.1333497 -13.09533868 3201.180
## 2 lcaT1Middle 2.05223339 0.2410168 8.51489685 2486.215
## 3 lcaT1High 3.31298684 0.2145886 15.43878319 3393.767
## 4 deprivationSIMD_2 0.01092666 0.1893618 0.05770256 3585.001
## 5 deprivationSIMD_3 0.07106083 0.1863199 0.38139143 3333.752
## 6 deprivationSIMD_4 -0.64834552 0.2149082 -3.01684933 3534.271
## 7 deprivationSIMD_5 -0.44423859 0.2078080 -2.13773583 2921.773
## 8 lcaT1Middle:deprivationSIMD_2 0.12866793 0.3408672 0.37747228 2340.421
## 9 lcaT1High:deprivationSIMD_2 -0.01633763 0.3082935 -0.05299377 3786.868
## 10 lcaT1Middle:deprivationSIMD_3 -0.12505026 0.3409031 -0.36682051 2661.730
## 11 lcaT1High:deprivationSIMD_3 0.04688613 0.3061281 0.15315854 3174.414
## 12 lcaT1Middle:deprivationSIMD_4 0.46230947 0.3470174 1.33223722 3316.506
## 13 lcaT1High:deprivationSIMD_4 0.43194974 0.3145081 1.37341372 3852.284
## 14 lcaT1Middle:deprivationSIMD_5 0.44409275 0.3478262 1.27676626 2524.584
## 15 lcaT1High:deprivationSIMD_5 0.35219871 0.3163132 1.11344922 2393.207
## p.value OR OR_lower OR_upper
## 1 3.301268e-38 0.1744251 0.1343072 0.2265263
## 2 2.840762e-17 7.7852692 4.8541811 12.4862290
## 3 5.065286e-52 27.4670428 18.0364226 41.8286075
## 4 9.539888e-01 1.0109866 0.6975211 1.4653232
## 5 7.029371e-01 1.0736465 0.7451824 1.5468923
## 6 2.572364e-03 0.5229102 0.3431576 0.7968207
## 7 3.262095e-02 0.6413124 0.4267563 0.9637388
## 8 7.058569e-01 1.1373124 0.5830794 2.2183590
## 9 9.577397e-01 0.9837951 0.5376256 1.8002356
## 10 7.137821e-01 0.8824526 0.4523856 1.7213689
## 11 8.782830e-01 1.0480027 0.5751497 1.9096064
## 12 1.828738e-01 1.5877366 0.8042504 3.1344807
## 13 1.697037e-01 1.5402577 0.8315315 2.8530415
## 14 2.018022e-01 1.5590751 0.7884813 3.0827809
## 15 2.656274e-01 1.4221911 0.7650797 2.6436820## Outcome: middle to high
mice.step3.mihiD <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "High" ~ lcaT1*deprivation, family = binomial,
subset = lcaT2 %in% c("Middle", "High")))
logreg.step3.mihiD <- summary(pool(mice.step3.mihiD))
logreg.step3.mihiD## term estimate std.error statistic df
## 1 (Intercept) -0.043414215 0.1739066 -0.249640985 1650.322
## 2 lcaT1Middle -0.039811289 0.2517969 -0.158108750 1523.182
## 3 lcaT1High 1.473634737 0.2383155 6.183544425 1402.372
## 4 deprivationSIMD_2 -0.238719563 0.2387201 -0.999997571 2196.206
## 5 deprivationSIMD_3 0.004999098 0.2429752 0.020574519 1612.702
## 6 deprivationSIMD_4 -0.102283901 0.2815389 -0.363302966 1999.808
## 7 deprivationSIMD_5 -0.293406777 0.2640270 -1.111275591 1702.267
## 8 lcaT1Middle:deprivationSIMD_2 0.172957441 0.3463979 0.499302835 1575.047
## 9 lcaT1High:deprivationSIMD_2 -0.054972655 0.3238750 -0.169734194 2228.954
## 10 lcaT1Middle:deprivationSIMD_3 0.106587672 0.3602919 0.295836987 1431.736
## 11 lcaT1High:deprivationSIMD_3 -0.161979661 0.3303424 -0.490338753 1231.301
## 12 lcaT1Middle:deprivationSIMD_4 -0.001425229 0.3816521 -0.003734366 1312.742
## 13 lcaT1High:deprivationSIMD_4 0.052506740 0.3626621 0.144781441 1916.319
## 14 lcaT1Middle:deprivationSIMD_5 0.375196692 0.3651187 1.027601840 1915.835
## 15 lcaT1High:deprivationSIMD_5 0.087989340 0.3490087 0.252112188 1069.884
## p.value
## 1 8.028961e-01
## 2 8.743921e-01
## 3 8.204792e-10
## 4 3.174218e-01
## 5 9.835876e-01
## 6 7.164170e-01
## 7 2.666067e-01
## 8 6.176358e-01
## 9 8.652346e-01
## 10 7.673975e-01
## 11 6.239816e-01
## 12 9.970210e-01
## 13 8.848987e-01
## 14 3.042669e-01
## 15 8.010028e-01coefs.mihiD <- logreg.step3.mihiD
coefs.mihiD$OR <- exp(coefs.mihiD$estimate)
coefs.mihiD$OR_lower <- exp(coefs.mihiD$estimate - 1.96*coefs.mihiD$std.error)
coefs.mihiD$OR_upper <- exp(coefs.mihiD$estimate + 1.96*coefs.mihiD$std.error)
coefs.mihiD## term estimate std.error statistic df
## 1 (Intercept) -0.043414215 0.1739066 -0.249640985 1650.322
## 2 lcaT1Middle -0.039811289 0.2517969 -0.158108750 1523.182
## 3 lcaT1High 1.473634737 0.2383155 6.183544425 1402.372
## 4 deprivationSIMD_2 -0.238719563 0.2387201 -0.999997571 2196.206
## 5 deprivationSIMD_3 0.004999098 0.2429752 0.020574519 1612.702
## 6 deprivationSIMD_4 -0.102283901 0.2815389 -0.363302966 1999.808
## 7 deprivationSIMD_5 -0.293406777 0.2640270 -1.111275591 1702.267
## 8 lcaT1Middle:deprivationSIMD_2 0.172957441 0.3463979 0.499302835 1575.047
## 9 lcaT1High:deprivationSIMD_2 -0.054972655 0.3238750 -0.169734194 2228.954
## 10 lcaT1Middle:deprivationSIMD_3 0.106587672 0.3602919 0.295836987 1431.736
## 11 lcaT1High:deprivationSIMD_3 -0.161979661 0.3303424 -0.490338753 1231.301
## 12 lcaT1Middle:deprivationSIMD_4 -0.001425229 0.3816521 -0.003734366 1312.742
## 13 lcaT1High:deprivationSIMD_4 0.052506740 0.3626621 0.144781441 1916.319
## 14 lcaT1Middle:deprivationSIMD_5 0.375196692 0.3651187 1.027601840 1915.835
## 15 lcaT1High:deprivationSIMD_5 0.087989340 0.3490087 0.252112188 1069.884
## p.value OR OR_lower OR_upper
## 1 8.028961e-01 0.9575147 0.6809468 1.346411
## 2 8.743921e-01 0.9609708 0.5866463 1.574142
## 3 8.204792e-10 4.3650722 2.7361075 6.963855
## 4 3.174218e-01 0.7876357 0.4933132 1.257558
## 5 9.835876e-01 1.0050116 0.6242325 1.618064
## 6 7.164170e-01 0.9027732 0.5199097 1.567579
## 7 2.666067e-01 0.7457187 0.4444580 1.251179
## 8 6.176358e-01 1.1888155 0.6029130 2.344090
## 9 8.652346e-01 0.9465110 0.5016927 1.785721
## 10 7.673975e-01 1.1124755 0.5490397 2.254120
## 11 6.239816e-01 0.8504585 0.4451025 1.624973
## 12 9.970210e-01 0.9985758 0.4726201 2.109842
## 13 8.848987e-01 1.0539097 0.5177251 2.145397
## 14 3.042669e-01 1.4552776 0.7114600 2.976742
## 15 8.010028e-01 1.0919765 0.5509740 2.164190## Outcome: middle to low
mice.step3.miloD <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Low" ~ lcaT1*deprivation, family = binomial,
subset = lcaT2 %in% c("Middle", "Low")))
logreg.step3.miloD <- summary(pool(mice.step3.miloD))
logreg.step3.miloD## term estimate std.error statistic df
## 1 (Intercept) 1.70284556 0.1311841 12.9805777 2607.5331
## 2 lcaT1Middle -2.09204468 0.2399319 -8.7193286 1375.5308
## 3 lcaT1High -1.83935211 0.2500233 -7.3567215 1936.5442
## 4 deprivationSIMD_2 -0.24964623 0.1778297 -1.4038498 2817.0596
## 5 deprivationSIMD_3 -0.06606173 0.1842109 -0.3586200 2393.2834
## 6 deprivationSIMD_4 0.54606173 0.2054523 2.6578519 2840.5023
## 7 deprivationSIMD_5 0.15083182 0.1891958 0.7972261 2827.2780
## 8 lcaT1Middle:deprivationSIMD_2 0.04428951 0.3330558 0.1329793 1323.7855
## 9 lcaT1High:deprivationSIMD_2 -0.03863502 0.3456650 -0.1117701 2613.8558
## 10 lcaT1Middle:deprivationSIMD_3 0.23163793 0.3470709 0.6674080 919.4983
## 11 lcaT1High:deprivationSIMD_3 -0.20886579 0.3487281 -0.5989360 2350.9783
## 12 lcaT1Middle:deprivationSIMD_4 -0.46373480 0.3410734 -1.3596336 1507.2538
## 13 lcaT1High:deprivationSIMD_4 -0.37944311 0.3579117 -1.0601584 2531.8854
## 14 lcaT1Middle:deprivationSIMD_5 -0.06889607 0.3358869 -0.2051169 2102.0969
## 15 lcaT1High:deprivationSIMD_5 -0.26420938 0.3541569 -0.7460236 1127.5467
## p.value
## 1 2.212678e-37
## 2 7.918348e-18
## 3 2.769699e-13
## 4 1.604738e-01
## 5 7.199110e-01
## 6 7.908113e-03
## 7 4.253866e-01
## 8 8.942300e-01
## 9 9.110143e-01
## 10 5.046790e-01
## 11 5.492733e-01
## 12 1.741493e-01
## 13 2.891737e-01
## 14 8.375007e-01
## 15 4.558087e-01coefs.miloD <- logreg.step3.miloD
coefs.miloD$OR <- exp(coefs.miloD$estimate)
coefs.miloD$OR_lower <- exp(coefs.miloD$estimate - 1.96*coefs.miloD$std.error)
coefs.miloD$OR_upper <- exp(coefs.miloD$estimate + 1.96*coefs.miloD$std.error)
coefs.miloD## term estimate std.error statistic df
## 1 (Intercept) 1.70284556 0.1311841 12.9805777 2607.5331
## 2 lcaT1Middle -2.09204468 0.2399319 -8.7193286 1375.5308
## 3 lcaT1High -1.83935211 0.2500233 -7.3567215 1936.5442
## 4 deprivationSIMD_2 -0.24964623 0.1778297 -1.4038498 2817.0596
## 5 deprivationSIMD_3 -0.06606173 0.1842109 -0.3586200 2393.2834
## 6 deprivationSIMD_4 0.54606173 0.2054523 2.6578519 2840.5023
## 7 deprivationSIMD_5 0.15083182 0.1891958 0.7972261 2827.2780
## 8 lcaT1Middle:deprivationSIMD_2 0.04428951 0.3330558 0.1329793 1323.7855
## 9 lcaT1High:deprivationSIMD_2 -0.03863502 0.3456650 -0.1117701 2613.8558
## 10 lcaT1Middle:deprivationSIMD_3 0.23163793 0.3470709 0.6674080 919.4983
## 11 lcaT1High:deprivationSIMD_3 -0.20886579 0.3487281 -0.5989360 2350.9783
## 12 lcaT1Middle:deprivationSIMD_4 -0.46373480 0.3410734 -1.3596336 1507.2538
## 13 lcaT1High:deprivationSIMD_4 -0.37944311 0.3579117 -1.0601584 2531.8854
## 14 lcaT1Middle:deprivationSIMD_5 -0.06889607 0.3358869 -0.2051169 2102.0969
## 15 lcaT1High:deprivationSIMD_5 -0.26420938 0.3541569 -0.7460236 1127.5467
## p.value OR OR_lower OR_upper
## 1 2.212678e-37 5.4895460 4.24492733 7.0990887
## 2 7.918348e-18 0.1234345 0.07712629 0.1975471
## 3 2.769699e-13 0.1589204 0.09735435 0.2594201
## 4 1.604738e-01 0.7790764 0.54980451 1.1039560
## 5 7.199110e-01 0.9360731 0.65238841 1.3431152
## 6 7.908113e-03 1.7264404 1.15416284 2.5824749
## 7 4.253866e-01 1.1628011 0.80252529 1.6848146
## 8 8.942300e-01 1.0452849 0.54416663 2.0078787
## 9 9.110143e-01 0.9621018 0.48863550 1.8943361
## 10 5.046790e-01 1.2606632 0.63850799 2.4890396
## 11 5.492733e-01 0.8115041 0.40968251 1.6074373
## 12 1.741493e-01 0.6289303 0.32231093 1.2272415
## 13 2.891737e-01 0.6842424 0.33927305 1.3799729
## 14 8.375007e-01 0.9334237 0.48324359 1.8029826
## 15 4.558087e-01 0.7678128 0.38352254 1.5371624## high to low
mice.step3.hiloD <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Low" ~ lcaT1*deprivation, family = binomial,
subset = lcaT2 %in% c("High", "Low")))
logreg.step3.hiloD <- summary(pool(mice.step3.hiloD))
logreg.step3.hiloD## term estimate std.error statistic df
## 1 (Intercept) 1.74625978 0.1333497 13.09533868 3201.180
## 2 lcaT1Middle -2.05223339 0.2410168 -8.51489685 2486.215
## 3 lcaT1High -3.31298684 0.2145886 -15.43878319 3393.767
## 4 deprivationSIMD_2 -0.01092666 0.1893618 -0.05770256 3585.001
## 5 deprivationSIMD_3 -0.07106083 0.1863199 -0.38139143 3333.752
## 6 deprivationSIMD_4 0.64834552 0.2149082 3.01684933 3534.271
## 7 deprivationSIMD_5 0.44423859 0.2078080 2.13773583 2921.773
## 8 lcaT1Middle:deprivationSIMD_2 -0.12866793 0.3408672 -0.37747228 2340.421
## 9 lcaT1High:deprivationSIMD_2 0.01633763 0.3082935 0.05299377 3786.868
## 10 lcaT1Middle:deprivationSIMD_3 0.12505026 0.3409031 0.36682051 2661.730
## 11 lcaT1High:deprivationSIMD_3 -0.04688613 0.3061281 -0.15315854 3174.414
## 12 lcaT1Middle:deprivationSIMD_4 -0.46230947 0.3470174 -1.33223722 3316.506
## 13 lcaT1High:deprivationSIMD_4 -0.43194974 0.3145081 -1.37341372 3852.284
## 14 lcaT1Middle:deprivationSIMD_5 -0.44409275 0.3478262 -1.27676626 2524.584
## 15 lcaT1High:deprivationSIMD_5 -0.35219871 0.3163132 -1.11344922 2393.207
## p.value
## 1 3.301268e-38
## 2 2.840762e-17
## 3 5.065286e-52
## 4 9.539888e-01
## 5 7.029371e-01
## 6 2.572364e-03
## 7 3.262095e-02
## 8 7.058569e-01
## 9 9.577397e-01
## 10 7.137821e-01
## 11 8.782830e-01
## 12 1.828738e-01
## 13 1.697037e-01
## 14 2.018022e-01
## 15 2.656274e-01coefs.hiloD <- logreg.step3.hiloD
coefs.hiloD$OR <- exp(coefs.hiloD$estimate)
coefs.hiloD$OR_lower <- exp(coefs.hiloD$estimate - 1.96*coefs.hiloD$std.error)
coefs.hiloD$OR_upper <- exp(coefs.hiloD$estimate + 1.96*coefs.hiloD$std.error)
coefs.hiloD## term estimate std.error statistic df
## 1 (Intercept) 1.74625978 0.1333497 13.09533868 3201.180
## 2 lcaT1Middle -2.05223339 0.2410168 -8.51489685 2486.215
## 3 lcaT1High -3.31298684 0.2145886 -15.43878319 3393.767
## 4 deprivationSIMD_2 -0.01092666 0.1893618 -0.05770256 3585.001
## 5 deprivationSIMD_3 -0.07106083 0.1863199 -0.38139143 3333.752
## 6 deprivationSIMD_4 0.64834552 0.2149082 3.01684933 3534.271
## 7 deprivationSIMD_5 0.44423859 0.2078080 2.13773583 2921.773
## 8 lcaT1Middle:deprivationSIMD_2 -0.12866793 0.3408672 -0.37747228 2340.421
## 9 lcaT1High:deprivationSIMD_2 0.01633763 0.3082935 0.05299377 3786.868
## 10 lcaT1Middle:deprivationSIMD_3 0.12505026 0.3409031 0.36682051 2661.730
## 11 lcaT1High:deprivationSIMD_3 -0.04688613 0.3061281 -0.15315854 3174.414
## 12 lcaT1Middle:deprivationSIMD_4 -0.46230947 0.3470174 -1.33223722 3316.506
## 13 lcaT1High:deprivationSIMD_4 -0.43194974 0.3145081 -1.37341372 3852.284
## 14 lcaT1Middle:deprivationSIMD_5 -0.44409275 0.3478262 -1.27676626 2524.584
## 15 lcaT1High:deprivationSIMD_5 -0.35219871 0.3163132 -1.11344922 2393.207
## p.value OR OR_lower OR_upper
## 1 3.301268e-38 5.73311937 4.41449886 7.44561472
## 2 2.840762e-17 0.12844771 0.08008823 0.20600797
## 3 5.065286e-52 0.03640727 0.02390708 0.05544337
## 4 9.539888e-01 0.98913282 0.68244328 1.43364841
## 5 7.029371e-01 0.93140523 0.64645743 1.34195334
## 6 2.572364e-03 1.91237423 1.25498741 2.91411307
## 7 3.262095e-02 1.55930248 1.03762550 2.34325797
## 8 7.058569e-01 0.87926589 0.45078367 1.71503220
## 9 9.577397e-01 1.01647182 0.55548286 1.86003034
## 10 7.137821e-01 1.13320540 0.58093299 2.21050363
## 11 8.782830e-01 0.95419605 0.52366813 1.73867770
## 12 1.828738e-01 0.62982740 0.31903211 1.24339382
## 13 1.697037e-01 0.64924201 0.35050314 1.20260031
## 14 2.018022e-01 0.64140593 0.32438244 1.26826090
## 15 2.656274e-01 0.70314038 0.37826032 1.30705329## high to mid
mice.step3.himiD <- with(data = ltadata.imputed,
exp = glm(lcaT2 == "Middle" ~ lcaT1*deprivation, family = binomial,
subset = lcaT2 %in% c("High", "Middle")))
logreg.step3.himiD <- summary(pool(mice.step3.himiD))
logreg.step3.himiD## term estimate std.error statistic df
## 1 (Intercept) 0.043414215 0.1739066 0.249640985 1650.322
## 2 lcaT1Middle 0.039811289 0.2517969 0.158108750 1523.182
## 3 lcaT1High -1.473634737 0.2383155 -6.183544425 1402.372
## 4 deprivationSIMD_2 0.238719563 0.2387201 0.999997571 2196.206
## 5 deprivationSIMD_3 -0.004999098 0.2429752 -0.020574519 1612.702
## 6 deprivationSIMD_4 0.102283901 0.2815389 0.363302966 1999.808
## 7 deprivationSIMD_5 0.293406777 0.2640270 1.111275591 1702.267
## 8 lcaT1Middle:deprivationSIMD_2 -0.172957441 0.3463979 -0.499302835 1575.047
## 9 lcaT1High:deprivationSIMD_2 0.054972655 0.3238750 0.169734194 2228.954
## 10 lcaT1Middle:deprivationSIMD_3 -0.106587672 0.3602919 -0.295836987 1431.736
## 11 lcaT1High:deprivationSIMD_3 0.161979661 0.3303424 0.490338753 1231.301
## 12 lcaT1Middle:deprivationSIMD_4 0.001425229 0.3816521 0.003734366 1312.742
## 13 lcaT1High:deprivationSIMD_4 -0.052506740 0.3626621 -0.144781441 1916.319
## 14 lcaT1Middle:deprivationSIMD_5 -0.375196692 0.3651187 -1.027601840 1915.835
## 15 lcaT1High:deprivationSIMD_5 -0.087989340 0.3490087 -0.252112188 1069.884
## p.value
## 1 8.028961e-01
## 2 8.743921e-01
## 3 8.204792e-10
## 4 3.174218e-01
## 5 9.835876e-01
## 6 7.164170e-01
## 7 2.666067e-01
## 8 6.176358e-01
## 9 8.652346e-01
## 10 7.673975e-01
## 11 6.239816e-01
## 12 9.970210e-01
## 13 8.848987e-01
## 14 3.042669e-01
## 15 8.010028e-01coefs.himiD <- logreg.step3.himiD
coefs.himiD$OR <- exp(coefs.himiD$estimate)
coefs.himiD$OR_lower <- exp(coefs.himiD$estimate - 1.96*coefs.himiD$std.error)
coefs.himiD$OR_upper <- exp(coefs.himiD$estimate + 1.96*coefs.himiD$std.error)
coefs.himiD## term estimate std.error statistic df
## 1 (Intercept) 0.043414215 0.1739066 0.249640985 1650.322
## 2 lcaT1Middle 0.039811289 0.2517969 0.158108750 1523.182
## 3 lcaT1High -1.473634737 0.2383155 -6.183544425 1402.372
## 4 deprivationSIMD_2 0.238719563 0.2387201 0.999997571 2196.206
## 5 deprivationSIMD_3 -0.004999098 0.2429752 -0.020574519 1612.702
## 6 deprivationSIMD_4 0.102283901 0.2815389 0.363302966 1999.808
## 7 deprivationSIMD_5 0.293406777 0.2640270 1.111275591 1702.267
## 8 lcaT1Middle:deprivationSIMD_2 -0.172957441 0.3463979 -0.499302835 1575.047
## 9 lcaT1High:deprivationSIMD_2 0.054972655 0.3238750 0.169734194 2228.954
## 10 lcaT1Middle:deprivationSIMD_3 -0.106587672 0.3602919 -0.295836987 1431.736
## 11 lcaT1High:deprivationSIMD_3 0.161979661 0.3303424 0.490338753 1231.301
## 12 lcaT1Middle:deprivationSIMD_4 0.001425229 0.3816521 0.003734366 1312.742
## 13 lcaT1High:deprivationSIMD_4 -0.052506740 0.3626621 -0.144781441 1916.319
## 14 lcaT1Middle:deprivationSIMD_5 -0.375196692 0.3651187 -1.027601840 1915.835
## 15 lcaT1High:deprivationSIMD_5 -0.087989340 0.3490087 -0.252112188 1069.884
## p.value OR OR_lower OR_upper
## 1 8.028961e-01 1.0443704 0.7427151 1.4685436
## 2 8.743921e-01 1.0406144 0.6352666 1.7046046
## 3 8.204792e-10 0.2290913 0.1435986 0.3654827
## 4 3.174218e-01 1.2696224 0.7951918 2.0271098
## 5 9.835876e-01 0.9950134 0.6180224 1.6019671
## 6 7.164170e-01 1.1076979 0.6379265 1.9234108
## 7 2.666067e-01 1.3409882 0.7992463 2.2499314
## 8 6.176358e-01 0.8411734 0.4266048 1.6586141
## 9 8.652346e-01 1.0565117 0.5599979 1.9932522
## 10 7.673975e-01 0.8988962 0.4436320 1.8213618
## 11 6.239816e-01 1.1758363 0.6153947 2.2466735
## 12 9.970210e-01 1.0014262 0.4739692 2.1158643
## 13 8.848987e-01 0.9488479 0.4661143 1.9315271
## 14 3.042669e-01 0.6871541 0.3359377 1.4055603
## 15 8.010028e-01 0.9157706 0.4620666 1.8149677Final Tables
Because we saved all our tables we can make a list of them and put them together in one dataframe, then turn that dataframe into a flextable in the same way we have done previously
#Make the list
all_tables <- list(
lomiG = coefs.lomiG, lohiG = coefs.lohiG,
mihiG = coefs.mihiG, miloG = coefs.miloG,
hiloG = coefs.hiloG, himiG = coefs.himiG,
lomiS = coefs.lomiS, lohiS = coefs.lohiS,
mihiS = coefs.mihiS, miloS = coefs.miloS,
hiloS = coefs.hiloS, himiS = coefs.himiS,
lomiD = coefs.lomiD, lohiD = coefs.lohiD,
mihiD = coefs.mihiD, miloD = coefs.miloD,
hiloD = coefs.hiloD, himiD = coefs.himiD
)
#Make the dataframe
combined_results <- bind_rows(
lapply(names(all_tables), function(x) {
df <- all_tables[[x]]
df$model <- x
df
})
)
combined_results <- combined_results %>%
dplyr::select(model, everything())
combined_results <- combined_results %>%
mutate(across(where(is.numeric), \(x) round(x, 3)))
#Make the flextable
ft <- flextable(combined_results)
ft <- autofit(ft)
ftmodel | term | estimate | std.error | statistic | df | p.value | OR | OR_lower | OR_upper |
|---|---|---|---|---|---|---|---|---|---|
lomiG | (Intercept) | -1.582 | 0.079 | -20.093 | 2,820.696 | 0.000 | 0.206 | 0.176 | 0.240 |
lomiG | lcaT1Middle | 2.093 | 0.146 | 14.311 | 2,535.052 | 0.000 | 8.107 | 6.087 | 10.799 |
lomiG | lcaT1High | 1.759 | 0.150 | 11.699 | 1,656.148 | 0.000 | 5.804 | 4.323 | 7.793 |
lomiG | genderMale | -0.379 | 0.121 | -3.120 | 1,987.209 | 0.002 | 0.685 | 0.540 | 0.869 |
lomiG | lcaT1Middle:genderMale | 0.096 | 0.210 | 0.454 | 2,459.888 | 0.650 | 1.100 | 0.729 | 1.661 |
lomiG | lcaT1High:genderMale | 0.510 | 0.221 | 2.309 | 2,099.554 | 0.021 | 1.665 | 1.080 | 2.567 |
lohiG | (Intercept) | -1.861 | 0.089 | -20.834 | 2,434.222 | 0.000 | 0.156 | 0.131 | 0.185 |
lohiG | lcaT1Middle | 2.260 | 0.157 | 14.435 | 1,518.171 | 0.000 | 9.582 | 7.050 | 13.023 |
lohiG | lcaT1High | 3.228 | 0.138 | 23.349 | 1,894.949 | 0.000 | 25.231 | 19.242 | 33.084 |
lohiG | genderMale | -0.121 | 0.129 | -0.943 | 2,353.689 | 0.346 | 0.886 | 0.688 | 1.140 |
lohiG | lcaT1Middle:genderMale | -0.104 | 0.219 | -0.473 | 1,755.154 | 0.636 | 0.901 | 0.586 | 1.386 |
lohiG | lcaT1High:genderMale | 0.439 | 0.200 | 2.190 | 2,038.980 | 0.029 | 1.550 | 1.047 | 2.295 |
mihiG | (Intercept) | -0.279 | 0.110 | -2.536 | 1,854.222 | 0.011 | 0.757 | 0.610 | 0.939 |
mihiG | lcaT1Middle | 0.167 | 0.156 | 1.069 | 1,597.759 | 0.285 | 1.182 | 0.870 | 1.606 |
mihiG | lcaT1High | 1.469 | 0.147 | 10.010 | 1,888.315 | 0.000 | 4.347 | 3.260 | 5.796 |
mihiG | genderMale | 0.258 | 0.165 | 1.560 | 1,515.505 | 0.119 | 1.294 | 0.936 | 1.789 |
mihiG | lcaT1Middle:genderMale | -0.199 | 0.229 | -0.871 | 1,801.388 | 0.384 | 0.819 | 0.523 | 1.283 |
mihiG | lcaT1High:genderMale | -0.071 | 0.217 | -0.329 | 1,309.721 | 0.742 | 0.931 | 0.609 | 1.424 |
miloG | (Intercept) | 1.582 | 0.079 | 20.093 | 2,820.696 | 0.000 | 4.865 | 4.169 | 5.677 |
miloG | lcaT1Middle | -2.093 | 0.146 | -14.311 | 2,535.052 | 0.000 | 0.123 | 0.093 | 0.164 |
miloG | lcaT1High | -1.759 | 0.150 | -11.699 | 1,656.148 | 0.000 | 0.172 | 0.128 | 0.231 |
miloG | genderMale | 0.379 | 0.121 | 3.120 | 1,987.209 | 0.002 | 1.461 | 1.151 | 1.854 |
miloG | lcaT1Middle:genderMale | -0.096 | 0.210 | -0.454 | 2,459.888 | 0.650 | 0.909 | 0.602 | 1.372 |
miloG | lcaT1High:genderMale | -0.510 | 0.221 | -2.309 | 2,099.554 | 0.021 | 0.601 | 0.390 | 0.926 |
hiloG | (Intercept) | 1.861 | 0.089 | 20.834 | 2,434.222 | 0.000 | 6.431 | 5.398 | 7.661 |
hiloG | lcaT1Middle | -2.260 | 0.157 | -14.435 | 1,518.171 | 0.000 | 0.104 | 0.077 | 0.142 |
hiloG | lcaT1High | -3.228 | 0.138 | -23.349 | 1,894.949 | 0.000 | 0.040 | 0.030 | 0.052 |
hiloG | genderMale | 0.121 | 0.129 | 0.943 | 2,353.689 | 0.346 | 1.129 | 0.877 | 1.453 |
hiloG | lcaT1Middle:genderMale | 0.104 | 0.219 | 0.473 | 1,755.154 | 0.636 | 1.109 | 0.722 | 1.706 |
hiloG | lcaT1High:genderMale | -0.439 | 0.200 | -2.190 | 2,038.980 | 0.029 | 0.645 | 0.436 | 0.955 |
himiG | (Intercept) | 0.279 | 0.110 | 2.536 | 1,854.222 | 0.011 | 1.322 | 1.065 | 1.640 |
himiG | lcaT1Middle | -0.167 | 0.156 | -1.069 | 1,597.759 | 0.285 | 0.846 | 0.623 | 1.150 |
himiG | lcaT1High | -1.469 | 0.147 | -10.010 | 1,888.315 | 0.000 | 0.230 | 0.173 | 0.307 |
himiG | genderMale | -0.258 | 0.165 | -1.560 | 1,515.505 | 0.119 | 0.773 | 0.559 | 1.068 |
himiG | lcaT1Middle:genderMale | 0.199 | 0.229 | 0.871 | 1,801.388 | 0.384 | 1.221 | 0.779 | 1.912 |
himiG | lcaT1High:genderMale | 0.071 | 0.217 | 0.329 | 1,309.721 | 0.742 | 1.074 | 0.702 | 1.642 |
lomiS | (Intercept) | -1.775 | 0.084 | -21.020 | 1,717.040 | 0.000 | 0.169 | 0.144 | 0.200 |
lomiS | lcaT1Middle | 2.154 | 0.148 | 14.536 | 2,401.199 | 0.000 | 8.619 | 6.447 | 11.524 |
lomiS | lcaT1High | 2.112 | 0.156 | 13.524 | 2,275.916 | 0.000 | 8.266 | 6.086 | 11.227 |
lomiS | scoreBeginner | 0.036 | 0.119 | 0.304 | 2,340.355 | 0.761 | 1.037 | 0.821 | 1.311 |
lomiS | lcaT1Middle:scoreBeginner | -0.059 | 0.208 | -0.285 | 2,699.549 | 0.776 | 0.942 | 0.626 | 1.418 |
lomiS | lcaT1High:scoreBeginner | -0.230 | 0.219 | -1.053 | 2,577.345 | 0.292 | 0.794 | 0.518 | 1.219 |
lohiS | (Intercept) | -1.877 | 0.087 | -21.472 | 2,721.902 | 0.000 | 0.153 | 0.129 | 0.182 |
lohiS | lcaT1Middle | 2.190 | 0.153 | 14.329 | 2,014.613 | 0.000 | 8.935 | 6.622 | 12.056 |
lohiS | lcaT1High | 3.387 | 0.141 | 24.077 | 3,235.221 | 0.000 | 29.566 | 22.442 | 38.951 |
lohiS | scoreBeginner | -0.094 | 0.128 | -0.735 | 3,155.345 | 0.462 | 0.910 | 0.708 | 1.170 |
lohiS | lcaT1Middle:scoreBeginner | 0.034 | 0.219 | 0.155 | 1,715.399 | 0.877 | 1.035 | 0.673 | 1.591 |
lohiS | lcaT1High:scoreBeginner | 0.120 | 0.197 | 0.608 | 3,522.629 | 0.543 | 1.128 | 0.766 | 1.660 |
mihiS | (Intercept) | -0.102 | 0.113 | -0.898 | 1,510.577 | 0.369 | 0.903 | 0.723 | 1.128 |
mihiS | lcaT1Middle | 0.036 | 0.160 | 0.225 | 1,503.794 | 0.822 | 1.037 | 0.758 | 1.419 |
mihiS | lcaT1High | 1.274 | 0.148 | 8.610 | 2,128.887 | 0.000 | 3.577 | 2.676 | 4.781 |
mihiS | scoreBeginner | -0.130 | 0.164 | -0.795 | 1,620.012 | 0.427 | 0.878 | 0.637 | 1.210 |
mihiS | lcaT1Middle:scoreBeginner | 0.093 | 0.231 | 0.405 | 1,217.697 | 0.686 | 1.098 | 0.698 | 1.727 |
mihiS | lcaT1High:scoreBeginner | 0.350 | 0.212 | 1.652 | 2,268.184 | 0.099 | 1.419 | 0.937 | 2.151 |
miloS | (Intercept) | 1.775 | 0.084 | 21.020 | 1,717.040 | 0.000 | 5.900 | 5.000 | 6.962 |
miloS | lcaT1Middle | -2.154 | 0.148 | -14.536 | 2,401.199 | 0.000 | 0.116 | 0.087 | 0.155 |
miloS | lcaT1High | -2.112 | 0.156 | -13.524 | 2,275.916 | 0.000 | 0.121 | 0.089 | 0.164 |
miloS | scoreBeginner | -0.036 | 0.119 | -0.304 | 2,340.355 | 0.761 | 0.964 | 0.763 | 1.219 |
miloS | lcaT1Middle:scoreBeginner | 0.059 | 0.208 | 0.285 | 2,699.549 | 0.776 | 1.061 | 0.705 | 1.597 |
miloS | lcaT1High:scoreBeginner | 0.230 | 0.219 | 1.053 | 2,577.345 | 0.292 | 1.259 | 0.820 | 1.932 |
hiloS | (Intercept) | 1.877 | 0.087 | 21.472 | 2,721.902 | 0.000 | 6.532 | 5.503 | 7.752 |
hiloS | lcaT1Middle | -2.190 | 0.153 | -14.329 | 2,014.613 | 0.000 | 0.112 | 0.083 | 0.151 |
hiloS | lcaT1High | -3.387 | 0.141 | -24.077 | 3,235.221 | 0.000 | 0.034 | 0.026 | 0.045 |
hiloS | scoreBeginner | 0.094 | 0.128 | 0.735 | 3,155.345 | 0.462 | 1.099 | 0.855 | 1.412 |
hiloS | lcaT1Middle:scoreBeginner | -0.034 | 0.219 | -0.155 | 1,715.399 | 0.877 | 0.966 | 0.629 | 1.486 |
hiloS | lcaT1High:scoreBeginner | -0.120 | 0.197 | -0.608 | 3,522.629 | 0.543 | 0.887 | 0.602 | 1.306 |
himiS | (Intercept) | 0.102 | 0.113 | 0.898 | 1,510.577 | 0.369 | 1.107 | 0.887 | 1.382 |
himiS | lcaT1Middle | -0.036 | 0.160 | -0.225 | 1,503.794 | 0.822 | 0.965 | 0.705 | 1.320 |
himiS | lcaT1High | -1.274 | 0.148 | -8.610 | 2,128.887 | 0.000 | 0.280 | 0.209 | 0.374 |
himiS | scoreBeginner | 0.130 | 0.164 | 0.795 | 1,620.012 | 0.427 | 1.139 | 0.826 | 1.571 |
himiS | lcaT1Middle:scoreBeginner | -0.093 | 0.231 | -0.405 | 1,217.697 | 0.686 | 0.911 | 0.579 | 1.432 |
himiS | lcaT1High:scoreBeginner | -0.350 | 0.212 | -1.652 | 2,268.184 | 0.099 | 0.704 | 0.465 | 1.068 |
lomiD | (Intercept) | -1.703 | 0.131 | -12.981 | 2,607.533 | 0.000 | 0.182 | 0.141 | 0.236 |
lomiD | lcaT1Middle | 2.092 | 0.240 | 8.719 | 1,375.531 | 0.000 | 8.101 | 5.062 | 12.966 |
lomiD | lcaT1High | 1.839 | 0.250 | 7.357 | 1,936.544 | 0.000 | 6.292 | 3.855 | 10.272 |
lomiD | deprivationSIMD_2 | 0.250 | 0.178 | 1.404 | 2,817.060 | 0.160 | 1.284 | 0.906 | 1.819 |
lomiD | deprivationSIMD_3 | 0.066 | 0.184 | 0.359 | 2,393.283 | 0.720 | 1.068 | 0.745 | 1.533 |
lomiD | deprivationSIMD_4 | -0.546 | 0.205 | -2.658 | 2,840.502 | 0.008 | 0.579 | 0.387 | 0.866 |
lomiD | deprivationSIMD_5 | -0.151 | 0.189 | -0.797 | 2,827.278 | 0.425 | 0.860 | 0.594 | 1.246 |
lomiD | lcaT1Middle:deprivationSIMD_2 | -0.044 | 0.333 | -0.133 | 1,323.786 | 0.894 | 0.957 | 0.498 | 1.838 |
lomiD | lcaT1High:deprivationSIMD_2 | 0.039 | 0.346 | 0.112 | 2,613.856 | 0.911 | 1.039 | 0.528 | 2.047 |
lomiD | lcaT1Middle:deprivationSIMD_3 | -0.232 | 0.347 | -0.667 | 919.498 | 0.505 | 0.793 | 0.402 | 1.566 |
lomiD | lcaT1High:deprivationSIMD_3 | 0.209 | 0.349 | 0.599 | 2,350.978 | 0.549 | 1.232 | 0.622 | 2.441 |
lomiD | lcaT1Middle:deprivationSIMD_4 | 0.464 | 0.341 | 1.360 | 1,507.254 | 0.174 | 1.590 | 0.815 | 3.103 |
lomiD | lcaT1High:deprivationSIMD_4 | 0.379 | 0.358 | 1.060 | 2,531.885 | 0.289 | 1.461 | 0.725 | 2.947 |
lomiD | lcaT1Middle:deprivationSIMD_5 | 0.069 | 0.336 | 0.205 | 2,102.097 | 0.838 | 1.071 | 0.555 | 2.069 |
lomiD | lcaT1High:deprivationSIMD_5 | 0.264 | 0.354 | 0.746 | 1,127.547 | 0.456 | 1.302 | 0.651 | 2.607 |
lohiD | (Intercept) | -1.746 | 0.133 | -13.095 | 3,201.180 | 0.000 | 0.174 | 0.134 | 0.227 |
lohiD | lcaT1Middle | 2.052 | 0.241 | 8.515 | 2,486.215 | 0.000 | 7.785 | 4.854 | 12.486 |
lohiD | lcaT1High | 3.313 | 0.215 | 15.439 | 3,393.767 | 0.000 | 27.467 | 18.036 | 41.829 |
lohiD | deprivationSIMD_2 | 0.011 | 0.189 | 0.058 | 3,585.001 | 0.954 | 1.011 | 0.698 | 1.465 |
lohiD | deprivationSIMD_3 | 0.071 | 0.186 | 0.381 | 3,333.752 | 0.703 | 1.074 | 0.745 | 1.547 |
lohiD | deprivationSIMD_4 | -0.648 | 0.215 | -3.017 | 3,534.271 | 0.003 | 0.523 | 0.343 | 0.797 |
lohiD | deprivationSIMD_5 | -0.444 | 0.208 | -2.138 | 2,921.773 | 0.033 | 0.641 | 0.427 | 0.964 |
lohiD | lcaT1Middle:deprivationSIMD_2 | 0.129 | 0.341 | 0.377 | 2,340.421 | 0.706 | 1.137 | 0.583 | 2.218 |
lohiD | lcaT1High:deprivationSIMD_2 | -0.016 | 0.308 | -0.053 | 3,786.868 | 0.958 | 0.984 | 0.538 | 1.800 |
lohiD | lcaT1Middle:deprivationSIMD_3 | -0.125 | 0.341 | -0.367 | 2,661.730 | 0.714 | 0.882 | 0.452 | 1.721 |
lohiD | lcaT1High:deprivationSIMD_3 | 0.047 | 0.306 | 0.153 | 3,174.414 | 0.878 | 1.048 | 0.575 | 1.910 |
lohiD | lcaT1Middle:deprivationSIMD_4 | 0.462 | 0.347 | 1.332 | 3,316.506 | 0.183 | 1.588 | 0.804 | 3.134 |
lohiD | lcaT1High:deprivationSIMD_4 | 0.432 | 0.315 | 1.373 | 3,852.284 | 0.170 | 1.540 | 0.832 | 2.853 |
lohiD | lcaT1Middle:deprivationSIMD_5 | 0.444 | 0.348 | 1.277 | 2,524.584 | 0.202 | 1.559 | 0.788 | 3.083 |
lohiD | lcaT1High:deprivationSIMD_5 | 0.352 | 0.316 | 1.113 | 2,393.207 | 0.266 | 1.422 | 0.765 | 2.644 |
mihiD | (Intercept) | -0.043 | 0.174 | -0.250 | 1,650.322 | 0.803 | 0.958 | 0.681 | 1.346 |
mihiD | lcaT1Middle | -0.040 | 0.252 | -0.158 | 1,523.182 | 0.874 | 0.961 | 0.587 | 1.574 |
mihiD | lcaT1High | 1.474 | 0.238 | 6.184 | 1,402.372 | 0.000 | 4.365 | 2.736 | 6.964 |
mihiD | deprivationSIMD_2 | -0.239 | 0.239 | -1.000 | 2,196.206 | 0.317 | 0.788 | 0.493 | 1.258 |
mihiD | deprivationSIMD_3 | 0.005 | 0.243 | 0.021 | 1,612.702 | 0.984 | 1.005 | 0.624 | 1.618 |
mihiD | deprivationSIMD_4 | -0.102 | 0.282 | -0.363 | 1,999.808 | 0.716 | 0.903 | 0.520 | 1.568 |
mihiD | deprivationSIMD_5 | -0.293 | 0.264 | -1.111 | 1,702.267 | 0.267 | 0.746 | 0.444 | 1.251 |
mihiD | lcaT1Middle:deprivationSIMD_2 | 0.173 | 0.346 | 0.499 | 1,575.047 | 0.618 | 1.189 | 0.603 | 2.344 |
mihiD | lcaT1High:deprivationSIMD_2 | -0.055 | 0.324 | -0.170 | 2,228.954 | 0.865 | 0.947 | 0.502 | 1.786 |
mihiD | lcaT1Middle:deprivationSIMD_3 | 0.107 | 0.360 | 0.296 | 1,431.736 | 0.767 | 1.112 | 0.549 | 2.254 |
mihiD | lcaT1High:deprivationSIMD_3 | -0.162 | 0.330 | -0.490 | 1,231.301 | 0.624 | 0.850 | 0.445 | 1.625 |
mihiD | lcaT1Middle:deprivationSIMD_4 | -0.001 | 0.382 | -0.004 | 1,312.742 | 0.997 | 0.999 | 0.473 | 2.110 |
mihiD | lcaT1High:deprivationSIMD_4 | 0.053 | 0.363 | 0.145 | 1,916.319 | 0.885 | 1.054 | 0.518 | 2.145 |
mihiD | lcaT1Middle:deprivationSIMD_5 | 0.375 | 0.365 | 1.028 | 1,915.835 | 0.304 | 1.455 | 0.711 | 2.977 |
mihiD | lcaT1High:deprivationSIMD_5 | 0.088 | 0.349 | 0.252 | 1,069.884 | 0.801 | 1.092 | 0.551 | 2.164 |
miloD | (Intercept) | 1.703 | 0.131 | 12.981 | 2,607.533 | 0.000 | 5.490 | 4.245 | 7.099 |
miloD | lcaT1Middle | -2.092 | 0.240 | -8.719 | 1,375.531 | 0.000 | 0.123 | 0.077 | 0.198 |
miloD | lcaT1High | -1.839 | 0.250 | -7.357 | 1,936.544 | 0.000 | 0.159 | 0.097 | 0.259 |
miloD | deprivationSIMD_2 | -0.250 | 0.178 | -1.404 | 2,817.060 | 0.160 | 0.779 | 0.550 | 1.104 |
miloD | deprivationSIMD_3 | -0.066 | 0.184 | -0.359 | 2,393.283 | 0.720 | 0.936 | 0.652 | 1.343 |
miloD | deprivationSIMD_4 | 0.546 | 0.205 | 2.658 | 2,840.502 | 0.008 | 1.726 | 1.154 | 2.582 |
miloD | deprivationSIMD_5 | 0.151 | 0.189 | 0.797 | 2,827.278 | 0.425 | 1.163 | 0.803 | 1.685 |
miloD | lcaT1Middle:deprivationSIMD_2 | 0.044 | 0.333 | 0.133 | 1,323.786 | 0.894 | 1.045 | 0.544 | 2.008 |
miloD | lcaT1High:deprivationSIMD_2 | -0.039 | 0.346 | -0.112 | 2,613.856 | 0.911 | 0.962 | 0.489 | 1.894 |
miloD | lcaT1Middle:deprivationSIMD_3 | 0.232 | 0.347 | 0.667 | 919.498 | 0.505 | 1.261 | 0.639 | 2.489 |
miloD | lcaT1High:deprivationSIMD_3 | -0.209 | 0.349 | -0.599 | 2,350.978 | 0.549 | 0.812 | 0.410 | 1.607 |
miloD | lcaT1Middle:deprivationSIMD_4 | -0.464 | 0.341 | -1.360 | 1,507.254 | 0.174 | 0.629 | 0.322 | 1.227 |
miloD | lcaT1High:deprivationSIMD_4 | -0.379 | 0.358 | -1.060 | 2,531.885 | 0.289 | 0.684 | 0.339 | 1.380 |
miloD | lcaT1Middle:deprivationSIMD_5 | -0.069 | 0.336 | -0.205 | 2,102.097 | 0.838 | 0.933 | 0.483 | 1.803 |
miloD | lcaT1High:deprivationSIMD_5 | -0.264 | 0.354 | -0.746 | 1,127.547 | 0.456 | 0.768 | 0.384 | 1.537 |
hiloD | (Intercept) | 1.746 | 0.133 | 13.095 | 3,201.180 | 0.000 | 5.733 | 4.414 | 7.446 |
hiloD | lcaT1Middle | -2.052 | 0.241 | -8.515 | 2,486.215 | 0.000 | 0.128 | 0.080 | 0.206 |
hiloD | lcaT1High | -3.313 | 0.215 | -15.439 | 3,393.767 | 0.000 | 0.036 | 0.024 | 0.055 |
hiloD | deprivationSIMD_2 | -0.011 | 0.189 | -0.058 | 3,585.001 | 0.954 | 0.989 | 0.682 | 1.434 |
hiloD | deprivationSIMD_3 | -0.071 | 0.186 | -0.381 | 3,333.752 | 0.703 | 0.931 | 0.646 | 1.342 |
hiloD | deprivationSIMD_4 | 0.648 | 0.215 | 3.017 | 3,534.271 | 0.003 | 1.912 | 1.255 | 2.914 |
hiloD | deprivationSIMD_5 | 0.444 | 0.208 | 2.138 | 2,921.773 | 0.033 | 1.559 | 1.038 | 2.343 |
hiloD | lcaT1Middle:deprivationSIMD_2 | -0.129 | 0.341 | -0.377 | 2,340.421 | 0.706 | 0.879 | 0.451 | 1.715 |
hiloD | lcaT1High:deprivationSIMD_2 | 0.016 | 0.308 | 0.053 | 3,786.868 | 0.958 | 1.016 | 0.555 | 1.860 |
hiloD | lcaT1Middle:deprivationSIMD_3 | 0.125 | 0.341 | 0.367 | 2,661.730 | 0.714 | 1.133 | 0.581 | 2.211 |
hiloD | lcaT1High:deprivationSIMD_3 | -0.047 | 0.306 | -0.153 | 3,174.414 | 0.878 | 0.954 | 0.524 | 1.739 |
hiloD | lcaT1Middle:deprivationSIMD_4 | -0.462 | 0.347 | -1.332 | 3,316.506 | 0.183 | 0.630 | 0.319 | 1.243 |
hiloD | lcaT1High:deprivationSIMD_4 | -0.432 | 0.315 | -1.373 | 3,852.284 | 0.170 | 0.649 | 0.351 | 1.203 |
hiloD | lcaT1Middle:deprivationSIMD_5 | -0.444 | 0.348 | -1.277 | 2,524.584 | 0.202 | 0.641 | 0.324 | 1.268 |
hiloD | lcaT1High:deprivationSIMD_5 | -0.352 | 0.316 | -1.113 | 2,393.207 | 0.266 | 0.703 | 0.378 | 1.307 |
himiD | (Intercept) | 0.043 | 0.174 | 0.250 | 1,650.322 | 0.803 | 1.044 | 0.743 | 1.469 |
himiD | lcaT1Middle | 0.040 | 0.252 | 0.158 | 1,523.182 | 0.874 | 1.041 | 0.635 | 1.705 |
himiD | lcaT1High | -1.474 | 0.238 | -6.184 | 1,402.372 | 0.000 | 0.229 | 0.144 | 0.365 |
himiD | deprivationSIMD_2 | 0.239 | 0.239 | 1.000 | 2,196.206 | 0.317 | 1.270 | 0.795 | 2.027 |
himiD | deprivationSIMD_3 | -0.005 | 0.243 | -0.021 | 1,612.702 | 0.984 | 0.995 | 0.618 | 1.602 |
himiD | deprivationSIMD_4 | 0.102 | 0.282 | 0.363 | 1,999.808 | 0.716 | 1.108 | 0.638 | 1.923 |
himiD | deprivationSIMD_5 | 0.293 | 0.264 | 1.111 | 1,702.267 | 0.267 | 1.341 | 0.799 | 2.250 |
himiD | lcaT1Middle:deprivationSIMD_2 | -0.173 | 0.346 | -0.499 | 1,575.047 | 0.618 | 0.841 | 0.427 | 1.659 |
himiD | lcaT1High:deprivationSIMD_2 | 0.055 | 0.324 | 0.170 | 2,228.954 | 0.865 | 1.057 | 0.560 | 1.993 |
himiD | lcaT1Middle:deprivationSIMD_3 | -0.107 | 0.360 | -0.296 | 1,431.736 | 0.767 | 0.899 | 0.444 | 1.821 |
himiD | lcaT1High:deprivationSIMD_3 | 0.162 | 0.330 | 0.490 | 1,231.301 | 0.624 | 1.176 | 0.615 | 2.247 |
himiD | lcaT1Middle:deprivationSIMD_4 | 0.001 | 0.382 | 0.004 | 1,312.742 | 0.997 | 1.001 | 0.474 | 2.116 |
himiD | lcaT1High:deprivationSIMD_4 | -0.053 | 0.363 | -0.145 | 1,916.319 | 0.885 | 0.949 | 0.466 | 1.932 |
himiD | lcaT1Middle:deprivationSIMD_5 | -0.375 | 0.365 | -1.028 | 1,915.835 | 0.304 | 0.687 | 0.336 | 1.406 |
himiD | lcaT1High:deprivationSIMD_5 | -0.088 | 0.349 | -0.252 | 1,069.884 | 0.801 | 0.916 | 0.462 | 1.815 |
Then we can save that to a word document.
doc <- read_docx()
doc <- body_add_par(doc, "Regression Results", style = "heading 1")
doc <- body_add_flextable(doc, ft)
print(doc, target = "Regression_Results.docx")All we need to do is clean this up now taking the lines that match the T1 and T2 values we were interested in and discarding the rest as we discussed earlier. I have not found reliable code that can do this, but it is likely possible to automate!
And that’s all there is to it! Well done for conducting latent transition analysis!