# 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.

```{r}
3*7
```

So 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 "\<-"

```{r}
output<-3*7
```

That is now saved to "output." We can check this by calling it.

```{r}
output
```

We can also now use that saved object in another calculation, here I want to add ten to the output:

```{r}
output+10
```

```{r message=FALSE, warning=FALSE}
#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.

```{r}
#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()
```

Hopefully 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?

[You can view the presentation slides here](https://calummacgillivray.github.io/assets/files/lta_presentation.pdf)

# 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.

```{r message=FALSE, warning=FALSE}
#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.

```{r message=FALSE, warning=FALSE}
#Load Packages
library(tidyverse)
library(poLCA)
library(mice)
library(readxl)
library(officer)
library(flextable)
```

# Importing and Cleaning the Data

## Data Importing

```{r}
url <- "https://calummacgillivray.github.io/assets/files/LTA_dataset_1.xlsx"
destfile <- "LTA_dataset_1.xlsx"
curl::curl_download(url, destfile)
dataset_T1 <- read_excel(destfile)

#If they were stored in your working directory you would do something like this:
#dataset_T1 <- read_excel("dataset_T1.xlsx")
```
And time point 2:
```{r}
url <- "https://calummacgillivray.github.io/assets/files/LTA_dataset_2.xlsx"
destfile <- "LTA_dataset_2.xlsx"
curl::curl_download(url, destfile)
dataset_T2 <- read_excel(destfile)
```

## Data Cleaning

```{r}
#Inspecting the first dataset
names(dataset_T1)
summary(dataset_T1)
```

Lets loop a table for each column of the dataset, we will start with column 2 because 1 is ID.

```{r}
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
```

What do you notice?

Have a look at Y1_T1 and Y2_T1.

Now lets run it for the second dataset.

```{r}
for(i in 2:ncol(dataset_T2)) {
  print(table(dataset_T2[, i], useNA = "ifany"))
  cat("\n")} 
```

This 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:

```{r}
head(dataset_T1)
```

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.

```{r}
#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.

```{r}
head(dataset_merged)
```

Now we are in a good place to clean this up in a way that stops us having to constantly look at our codebook!

```{r}
#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.

```{r}
for(i in 2:ncol(dataset_merged)) {
  print(table(dataset_merged[, i], useNA = "ifany"))
  cat("\n")}
```

We 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.

```{r}
#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)
```

Now let's flip the positive variables as discussed.

```{r}
#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)](https://doi.org/10.1371/journal.pone.0317617) 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](https://learningstatisticswithr.com/book/descriptives.html).

I also highly recommend checking out this resource on automating your tables by [Leo Westbury](https://repository.ncrm.ac.uk/resources/online/all/producing_automated_outputs_r), 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.

```{r}
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 for
```

We should also set a seed to allow our analyses based on random number generation to be reproducible.

```{r}
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.

```{r}
#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)
```

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.

```{r}
lca_t1.3 <- poLCA(Indicators_T1, dataset_merged, nclass = 3, nrep = 10)
```

Again no error message, so let's keep going.

```{r}
lca_t1.4 <- poLCA(Indicators_T1, dataset_merged, nclass = 4, nrep = 10) # 4 classes
```

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.

```{r}
#We use maxiter to increase the maximum iterations
lca_t1.4 <- poLCA(Indicators_T1, dataset_merged, nclass = 4, nrep = 10, maxiter = 5000) # 4 classes

#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 classes
```

### Try 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.

```{r}
#Here are the rest for Time point 2 as well
lca_t2.2 <- poLCA(, dataset_merged, nclass = , nrep = 10) # 2 classes
lca_t2.3 <- poLCA(, dataset_merged, nclass = , nrep = 10, maxiter = 3000) # 3 classes
lca_t2.4 <- poLCA(, dataset_merged, nclass = , nrep = 10, maxiter = 6000) # 4 classes
#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 classes
```

To compare the fit of the models we will use a helper function that can harvest the relevant information from the outputs we stored.

```{r}
#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.

```{r}
#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)
#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.

```{r}
lc_t2.IC <- t(as.data.frame(lapply(list(, , ), ICfct)))
row.names(lc_t2.IC) <- 2:4
round(lc_t2.IC, 2)

lc_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).

```{r}
#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.

```{r}
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.

```{r}
#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.

```{r}
#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 table
```

And again for time point 2

```{r}
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_t2
```

Finally we will use the officer package to export these to a single word document.

```{r}
#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.

```{r}
#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

```{r}
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
```

Now again for T2

```{r}
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
```

Now one more dataset to make for a combined graph

```{r}
#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.

```{r fig.width=8}
#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.

```{r fig.width=8}
#Then we do the same from T2
probplot_T2 <- ggplot(data=, 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
```
Then we can see what they look like overlayed.
```{r fig.width=8}
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.

```{r}
high_probs_1 <- sapply(lca_t1.3$probs, function(item) item[, 2])
class_score_1 <- rowMeans(high_probs_1)
class_score_1

high_probs_2 <- sapply(lca_t2.3$probs, function(item) item[, 2])
class_score_2 <- rowMeans(high_probs_2)
class_score_2
```

We 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.

```{r}
#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)
table(lcaT2)
```

Now lets make the tables.

```{r}
#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_T1
```

And for T2.

```{r}
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_T2
```

Then like before, save them together as a word document.

```{r}
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](https://epirhandbook.com/en/new_pages/missing_data.html#assess-missingness-in-a-data-frame). 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.

```{r}
#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
```

```{r}
#Just a quick check
str(ltadata.sim)
```

Let's make deprivation a factor with levels like gender and score.

```{r}
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)
```

Now we can make the imputed datasets and save them to ltadata.imputed.

```{r}
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.
```

## Calculating 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.

```{r}
#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
```

We 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.

```{r}
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
```

We 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.

```{r}
ltadata.mice_mid <- with(data = ltadata.imputed, exp = glm(lcaT2 == "" ~ lcaT1 - 1, family = binomial))
logreg_mid <- summary(pool(ltadata.mice_mid))
logreg_mid          

coef_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
```

### Low class as outcome

### Try it yourself! #5

And do the same here but for "Low".

```{r}
ltadata.mice_low <- with(data = ltadata.imputed, exp = glm(lcaT2 == "" ~ lcaT1 - 1, family = binomial))
logreg_low <- summary(pool(ltadata.mice_low))
logreg_low

coef_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
```

## Make the Table

To interpret our step 2 findings, lets bring them together. We will be making another flextable

```{r}
#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 look
```

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.

```{r}
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.

```{r}
## 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.

```{r}
#Pool the results
logreg.step3.lomiG <- summary(pool(mice.step3.lomiG))
logreg.step3.lomiG

#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
```

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). 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.

```{r}
## Outcome: low to high
mice.step3.lohiG <- with(data = ltadata.imputed,
                         exp = glm(lcaT2 == "" ~ lcaT1*gender, family = binomial,
                                   subset = lcaT2 %in% c("", "") & gender %in% c("Male", "Female")))
logreg.step3.lohiG <- summary(pool(mice.step3.lohiG))
logreg.step3.lohiG

coefs.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
```

## Gender

Here is the full code for every transition.

```{r}
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

coefs.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

## 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

coefs.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

## 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

coefs.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

## 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

coefs.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

## 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

coefs.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

## 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

coefs.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
```

## Score

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.

```{r}
## 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

coefs.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

## 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

coefs.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

## 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

coefs.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

## 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

coefs.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

## 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

coefs.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

## 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

coefs.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
```

## Deprivation 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.

```{r}
## 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

coefs.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

## 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

coefs.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

## 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

coefs.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

## 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

coefs.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

## 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

coefs.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

## 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

coefs.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
```

## Final 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

```{r}
#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)
ft
```

Then we can save that to a word document.

```{r}
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!
