---
title: "SEM with R lavaan"
output: html_document
date: 2023.03.22
---

```{r setup, include=FALSE, message = FALSE, warning = FALSE}
knitr::opts_chunk$set(echo = TRUE)

# set custom functions for displaying 0-3 decimals in inline text

dec0 <- function(x){
        result <- paste(sprintf("%1.0f", x), collapse = ", ")
        return(result)
}

dec1 <- function(x){
        result <- paste(sprintf("%1.1f", x), collapse = ", ")
        return(result)
}

dec2 <- function(x){
        result <- paste(sprintf("%1.2f", x), collapse = ", ")
        return(result)
}

dec3 <- function(x){
        result <- paste(sprintf("%1.3f", x), collapse = ", ")
        return(result)
}
```

```{r install packages full set, message = FALSE, warning = FALSE, echo = FALSE}

# install.packages("compareGroups")
# install.packages("data.table")
# install.packages("devtools")
# install.packages("doBy")
# install.packages("extraoperators")
# install.packages("finalfit")
# install.packages("gapminder")
# install.packages("ggplot2")
# install.packages("ggpubr")
# install.packages("ggraph")
# install.packages("gmodels")
# install.packages("HLMdiag")
# install.packages("intervals")
# install.packages("JWileymisc")
# install.packages("lattice")
# install.packages("lme4")
# install.packages("lmerTest")
# install.packages("lmtest")
# install.packages("MatchIt")
# install.packages("mice")
# install.packages("mlmRev")
# install.packages("moments")
# install.packages("multilevelTools")
# install.packages("naniar")
# install.packages("nlme")
# install.packages("nullabor")
# install.packages("optmatch")
# install.packages("pan")
# install.packages("plyr")
# install.packages("psych")
# install.packages("RColorBrewer")
# install.packages("sandwich")
# install.packages("sjstats")
# install.packages("stargazer")
# install.packages("tableone")
# install.packages("texreg")
# install.packages("tidyr")
# install.packages("tidyverse")
# install.packages("tsibble")
# install.packages("apaTables")
# install.packages("sjPlot")
# install.packages("gridExtra")
# install.packages("emmeans")
# tinytex::install_tinytex() # knitting to pdf
# remotes::install_github("ha-pu/supportR")
# install.packages("simr")
```

```{r libraries full set, message = FALSE, warning = FALSE, echo = FALSE}

library(apaTables) # https://cran.r-project.org/web/packages/apaTables/apaTables.pdf
library(brms)
library(compareGroups)
library(compareGroups) # https://cran.r-project.org/web/packages/compareGroups/vignettes/compareGroups_vignette.html
library(data.table)
library(devtools)
library(doBy)
library(emmeans)
library(equatiomatic)
library(extraoperators)
library(finalfit)
library(gapminder)
library(ggmcmc)
library(ggplot2)
library(ggpubr)
library(ggraph)
library(ggridges)
library(ggthemes)
library(gmodels)
library(gridExtra)
library(HLMdiag)
library(intervals)
# library(JointAI)
library(JWileymisc)
library(lattice)
library(lme4)
library(lmerTest)
library(lmtest)
library(lubridate)
library(MatchIt)
# library(mice)
library(moments)
library(mlmRev)
library(multilevelTools)
library(naniar)
library(nlme)
library(nullabor)
library(optmatch)
library(pan)
library(parallel)
library(plyr)
library(psych)
library(RColorBrewer)
# library(rjags)
# library(rstan)
library(sandwich)
library(simr)
library(sjPlot)
library(sjstats)
library(stargazer)
library(stats)
library(supportR)
library(tableone)
library(texreg)
library(tidyr)
library(tidyverse)
library(tinytex)
library(tsibble)

```

```{r install packages SEM, message = FALSE}

# install.packages("lavaan", dependencies = TRUE)
# install.packages("qgraph", dependencies = TRUE)
# install.packages("semPlot", dependencies = TRUE)

```

```{r libraries SEM, message = FALSE}

library(lavaan)
library(qgraph)
library(semPlot)

```

  
## 1. Introduction   
  
[The lavaan package](https://lavaan.ugent.be/) is a versatile tool for Structural Equation Modeling (SEM) (or latent variable analysis, covariance structure analysis). SEM is a good (?) friend of a researcher who aims 1) to understand the patterns of correlation/covariance amongst the variables, and 2) to explain their variance with the model specified (see [Suhr, 2006](https://www.lexjansen.com/wuss/2006/tutorials/TUT-Suhr.pdf)).   

We will first look into regression and path models using only observed variables. Second, we will build a SEM that contains also latent (unobserved) variables.  
  
First (regression and path) models will be run with the familiar "dfttu" dataset that contains *autonomy*, *competence* and *score* variables alongside with "gender" and "group" variables. The SEM model will be run with another dataset which will be described later.  
  
If you use lavaan in your own analyses, please use the following reference:  
Yves Rosseel (2012). lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48(2), 1-36. URL http://www.jstatsoft.org/v48/i02/   
  
  
### 1.1. SEM terminology  
  
[Greek symbols](https://web.mit.edu/jmorzins/www/greek-alphabet.html)  
  
**observed/manifest variable** (something that you have measured: one column in the data frame, square box in diagram)  
  
**latent variable** (factor that is composed of more than one observed variables: does not exist in the data frame, circle in diagram)  
  
**exogenous variable** (independent: sends arrows, explains endogenous variable(s), may be observed (x) or latent ($\xi$)  
  
**endogenous variable** (dependent: receives arrows, is explained by other vars, may be observed (y) or latent ($\eta$)  


### 1.2 Basic model specification operators in lavaan  
  
"=~"   latent variable ("is measured by", e.g., SelfEfficacy =~ se1 + se2 + se3)   
  
"~"  regression ("is regressed on", e.g., score ~ autonomy + competence + gender)  
  
"~~"  (co)variance ("is correlated with", e.g., autonomy ~~ autonomy (variance) or autonomy ~~ competence (covariance))   
  
"~1"  intercept/mean (triangle in diagrams, e.g., autonomy ~ 1 estimates the mean of variable autonomy)  
  
"1*"   fixes parameter or loading to one (e.g., SelfEfficacy =~ 1*se1 + se2 + se3)  
  
NA*  frees parameter or loading (overrides default, e.g., SelfEfficacy =~ NA*se1 + se2 + se3)  
  
  
### 1.3 Estimators  
  
The *R lavaan* offers several estimator [options](https://lavaan.ugent.be/tutorial/est.html) (see also the  [manual](https://rdrr.io/cran/lavaan/man/lavOptions.html)), but we use in these examples one of the most common one, Maximum Likelihood (ML) estimation. As one of the names for "SEM" is "covariance structure modeling", it indicates that the task of ML (for details, see [Crisci, 2012](http://jaqm.ro/issues/volume-7,issue-2/pdfs/1_crisci.pdf)) is to compute the implied covariance matrix (specified in the model) and compare that against the actual covariance matrix based on the empirical data.  
  
  
  
## 2 Before the analysis  
  
Next we go through the steps that are needed before the estimation of a regression model.  
  
  
### 2.1 Set the working directory  
  
```{r set working directory, message = FALSE, warning = FALSE}

# where this current script is located?
current_dir <- dirname(rstudioapi::getSourceEditorContext()$path)

# set working directory to current script location
setwd(current_dir)

# check out the result
getwd()

```
  
  
### 2.2 Read and inspect the data  
  
```{r read the data and display summary plus six first rows, message = FALSE, warning = FALSE}

# read the data into data frame (df) called "dfttu"
dfttu <- read.table("dfttu_sem.txt", header=TRUE, sep= "\t")

# display summary of df
summary(dfttu)

# View(dfttu)

# display first six rows of df
head(dfttu)

```
  
```{r display descriptives stargazer, message = FALSE}

# install.packages("stargazer")
library(stargazer)
stargazer(dfttu, type = "text")

```

```{r display descriptives desctools, message = FALSE}

# install.packages("DescTools")
library(DescTools)
Desc(dfttu, plotit = TRUE)

```
  

```{r display descriptives psych, message = FALSE}

# install.packages("psych")
library(psych)
describe(dfttu)

```

  
The data (`r dec0(nrow(dfttu))` participants, i.e., rows) contains `r dec0(ncol(dfttu))` variables:  
  
group           0=control, 1=intervention  
gender          0=female, 1=male  
autonomy        1=totally disagree, ..., 5=totally agree (average of four items)  
competence      1=totally disagree, ..., 5=totally agree (average of four items)    
deep            1=totally disagree, ..., 5=totally agree (average of four items)    
organised       1=totally disagree, ..., 5=totally agree (average of four items)  
score           0-100 (mathematics score)  
  
We can see from the summary of the "dfttu" data frame (df) that it contains *no* missing values (NA). Usually we are happy about that, but in the real life NA's quite often do exist.  
  
Next we examine the distributions of the three central variables that we are going to use in the following exercises: Competence, autonomy and score.  
  
Frequency histograms show the count of each variable value:  
  
```{r histograms, message = FALSE}

hist(dfttu$competence, main="Histogram of Competence", xlab="Competence")
hist(dfttu$autonomy, main="Histogram of Autonomy", xlab="Autonomy")
hist(dfttu$score, main="Histogram of Score", xlab="Score")

```
  
  
Density plots of these three variables are based on kernel smoothing:  
    
  
```{r density plot competence autonomy score, message = FALSE, warning = FALSE}

# # basic plot
# plot(density(dfttu$competence))
# plot(density(dfttu$autonomy))
# plot(density(dfttu$score))

# "bells and whistles" version
# https://datavizpyr.com/add-vertical-line-to-density-plot-with-ggplot2/

# competence
# shorter way
mean_competence <- as.numeric(dec2(mean(dfttu$competence)))
median_competence <- as.numeric(dec2(median(dfttu$competence)))
mode_competence <- as.numeric(dec2(Mode(dfttu$competence))) # Mode from DescTools library
variance_competence <- as.numeric(dec2(var(dfttu$competence)))
sd_competence <- as.numeric(dec2(sd(dfttu$competence)))
skewness_competence <- as.numeric(dec2(skewness(dfttu$competence)))
kurtosis_competence <- as.numeric(dec2(kurtosis(dfttu$competence)))

# # longer way
# mean_competence <- dfttu %>% 
#   pull(competence) %>% 
#   mean() %>%
#   signif(3)

# display plot
ggplot(dfttu, aes(x=competence)) +
  geom_density( fill="dodgerblue", alpha=0.2) + 
  geom_vline(xintercept=mean_competence, size=1, color="red") +
  geom_text(aes(x=4.5, label=paste0("Mean (red)\n",mean_competence), y=1.5)) +
  geom_vline(xintercept=median_competence, size=1, color="blue") +
  geom_text(aes(x=4.5, label=paste0("Median (blue)\n",median_competence), y=1.2)) + 
  geom_vline(xintercept=mode_competence, size=1, color="orange", linetype="dashed") +
  geom_text(aes(x=4.5, label=paste0("Mode (orange)\n",mode_competence), y=0.9)) + 
  geom_text(aes(x=3, label=paste0("Variance\n",variance_competence), y=1.5)) + 
  geom_text(aes(x=3, label=paste0("Std dev\n",sd_competence), y=1.2)) + 
  geom_text(aes(x=3, label=paste0("Skewness\n",skewness_competence), y=0.9)) + 
  geom_text(aes(x=3, label=paste0("Kurtosis\n",kurtosis_competence), y=0.6))

# autonomy
mean_autonomy <- as.numeric(dec2(mean(dfttu$autonomy)))
median_autonomy <- as.numeric(dec2(median(dfttu$autonomy)))
mode_autonomy <- as.numeric(dec2(Mode(dfttu$autonomy))) # Mode from DescTools library
variance_autonomy <- as.numeric(dec2(var(dfttu$autonomy)))
sd_autonomy <- as.numeric(dec2(sd(dfttu$autonomy)))
skewness_autonomy <- as.numeric(dec2(skewness(dfttu$autonomy)))
kurtosis_autonomy <- as.numeric(dec2(kurtosis(dfttu$autonomy)))

ggplot(dfttu, aes(x=autonomy)) +
  geom_density( fill="dodgerblue", alpha=0.2) + 
  geom_vline(xintercept=mean_autonomy, size=1, color="red") +
  geom_text(aes(x=4.5, label=paste0("Mean (red)\n",mean_autonomy), y=1.5)) +
  geom_vline(xintercept=median_autonomy, size=1, color="blue") +
  geom_text(aes(x=4.5, label=paste0("Median (blue)\n",median_autonomy), y=1.2)) + 
  geom_vline(xintercept=mode_autonomy, size=1, color="orange", linetype="dashed") +
  geom_text(aes(x=4.5, label=paste0("Mode (orange)\n",mode_autonomy), y=0.9)) + 
  geom_text(aes(x=3, label=paste0("Variance\n",variance_autonomy), y=1.5)) + 
  geom_text(aes(x=3, label=paste0("Std dev\n",sd_autonomy), y=1.2)) + 
  geom_text(aes(x=3, label=paste0("Skewness\n",skewness_autonomy), y=0.9)) + 
  geom_text(aes(x=3, label=paste0("Kurtosis\n",kurtosis_autonomy), y=0.6))

# score
mean_score <- as.numeric(dec2(mean(dfttu$score)))
median_score <- as.numeric(dec2(median(dfttu$score)))
mode_score <- as.numeric(dec2(Mode(dfttu$score))) # Mode from DescTools library
variance_score <- as.numeric(dec2(var(dfttu$score)))
sd_score <- as.numeric(dec2(sd(dfttu$score)))
skewness_score <- as.numeric(dec2(skewness(dfttu$score)))
kurtosis_score <- as.numeric(dec2(kurtosis(dfttu$score)))

ggplot(dfttu, aes(x=score)) +
  geom_density( fill="dodgerblue", alpha=0.2) + 
  geom_vline(xintercept=mean_score, size=1, color="red") +
  geom_text(aes(x=75, label=paste0("Mean (red)\n",mean_score), y=0.01)) +
  geom_vline(xintercept=median_score, size=1, color="blue") +
  geom_text(aes(x=75, label=paste0("Median (blue)\n",median_score), y=0.007)) + 
  geom_vline(xintercept=mode_score, size=1, color="orange", linetype="dashed") +
  geom_text(aes(x=75, label=paste0("Mode (orange)\n",mode_score), y=0.004)) + 
  geom_text(aes(x=25, label=paste0("Variance\n",variance_score), y=0.01)) +
  geom_text(aes(x=25, label=paste0("Std dev\n",sd_score), y=0.007)) +
  geom_text(aes(x=25, label=paste0("Skewness\n",skewness_score), y=0.004)) +
  geom_text(aes(x=25, label=paste0("Kurtosis\n",kurtosis_score), y=0.001))

# install.packages("gridExtra")
# library(gridExtra)
# grid.arrange(plot_competence,plot_autonomy,plot_score, ncol=3)

```
  
Last figure of score density looks odd by the value of mode (11.11). This makes sense when we look at the frequencies of the original and rounded (zero digits) math score values: There are 17 occurences of 11.11 value.  
  
```{r math score frequencies, message = FALSE, warning = FALSE}

# table(round(dfttu$score, digits=0))
table(dfttu$score)

```
  
  
```{r kurtosis, message = FALSE, warning = FALSE}

# # install.packages("moments")
# library(moments)
# 
# # is the data more or less "tailed" than normal distribution ("bell curve")?
# kurtosis(dfttu)

```
  
  
  group      gender     autonomy   competence deep       organised  score  
  1.003331   1.736070   3.059942   2.721777   4.460560   2.844617   1.866836  
  
  
Kurtosis value of 3 indicates normal distribution (and it is called *mesokurtic*). If the value is less than 3, then the distribution is said to be *platykurtic*: It has thinner tails than a normal distribution, resulting in fewer extreme positive or negative values. If the value is greater than 3, then it is said to be *leptokurtic*: The shape is wider or flatter and the tails are fatter resulting in a greater chance of extreme positive or negative values. Competence and autonomy variables are quite close to (more or less) *normal* value of 3, but the score is clearly platycurtic.  
  

Skewness is a distribution shape measure that tells if the response values are concentrated more on the left (too many smaller values) or right (too many greater values) side of the mean.  
  
```{r skewness, message = FALSE, warning = FALSE}

# library(moments)
# 
# # is the data more towards low or high response values than the normal distribution?
# skewness(dfttu)

```
  
  
group       gender      autonomy    competence  deep        organised   score 
-0.05771633 -0.85794525 -0.17889672 -0.23497594 -0.65758103 -0.31181557 -0.10997036
  
  
Usually values greater than -1 or 1 indicate highly skewed (non-normal) data. Results above show that this is not the case here. Moderate skewness is attached to values -1 - -0.5 and 0.5 - 1. However, in our data autonomy, competence and score all have negative skewness values. It means that greater amount of values are larger than the mean value. Following density plots visualise this.  
    


  
```{r produce sqrt and log transformations on score, message = FALSE, warning = FALSE}

# sqrt
dfttu$scoreSqrt <- sqrt(dfttu$score)
# inverse sqrt
dfttu$scoreSqrtInv <- sqrt(max(dfttu$score+1)-dfttu$score)
# log10
dfttu$scoreLog <- log(dfttu$score)
# inverse log10
dfttu$scoreLogInv <- log(max(dfttu$score+1)-dfttu$score)

```
  
```{r sqrt and log transformations on score skewness, message = FALSE, warning = FALSE}

skewness(dfttu$score)
skewness(dfttu$scoreSqrt)
skewness(dfttu$scoreSqrtInv)
skewness(dfttu$scoreLog)
skewness(dfttu$scoreLogInv)

plot(density(dfttu$score))
plot(density(dfttu$scoreSqrt))
plot(density(dfttu$scoreSqrtInv))
plot(density(dfttu$scoreLog))
plot(density(dfttu$scoreLogInv))

```
  
Skewness is not getting better by transforming the score (as skewness values go further from the desired zero value). Density plots below show that the square root transformation is milder than logarithmic transformation.  
  
```{r sqrt and log transformations on score kurtosis, message = FALSE, warning = FALSE}

kurtosis(dfttu$score)
kurtosis(dfttu$scoreSqrt)
kurtosis(dfttu$scoreSqrtInv)
kurtosis(dfttu$scoreLog)
kurtosis(dfttu$scoreLogInv)

plot(density(dfttu$score))
plot(density(dfttu$scoreSqrt))
plot(density(dfttu$scoreSqrtInv))
plot(density(dfttu$scoreLog))
plot(density(dfttu$scoreLogInv))

```

  
The kurtosis of score benefits most from the square root transformation as the kurtosis value gets close to 3 (2.92, second figure). However, we decide to go on without transforming the score variable.    
  
  
### 2.3 Deal with the missing data  
  
We create a second df ("dfttuNA") and add approximately 40% of NA's in the *autonomy* and *organised* variables (in that new df).  
  
```{r add some missing values into new df, message = FALSE, warning = FALSE}

# create a copy of the "whole"no missing data" df ("dfttu")
dfttuNA <- dfttu

# create a new temporary df that only has autonomy and organised columns
# one way is to use "dplyr"
# library(dplyr)
# df_temp <- select(dfttuNA, autonomy, organised)

# but we opt to use "base R"
df_temp <- dfttu[c("autonomy","competence")]

# add appr. 40% of NA's into these two columns
set.seed(2021) # reproduce results
# https://stackoverflow.com/questions/27454265/randomly-insert-nas-into-dataframe-proportionaly/27454361#27454361
df_temp <- as.data.frame(lapply(df_temp,
                     function(cc) cc[ sample(c(TRUE, NA), prob = c(0.60, 0.40), size = length(cc), replace = TRUE) ]
                     )
              )

# overwrite dfttuNA autonomy and organised columns with the same vars containing NA's 
dfttuNA$autonomy <- as.numeric(df_temp[,1]) # first column of df_temp is autonomy
dfttuNA$organised <- as.numeric(df_temp[,2]) # seconf column of df_temp is organised

# display summary of the new df with NA's
summary(dfttuNA)
head(dfttuNA)

# write new df
write.table (dfttuNA, file = "dfttu_sem_na.txt", row.names = FALSE, sep= "\t")

```
  
Here we have it - 79 missing values in autonomy and 82 in organised - great! Let's observe the data patterns:  
  
```{r show missing data patterns, message = FALSE, warning = FALSE}

library(mice)
md.pattern(dfttuNA, rotate.names = TRUE)

```

  
First pattern is the complete data (75 participants) while the second one is the missing data in *organised* variable (54 participants). The third pattern represents the missigness in *autonomy* (51 participants), and finally, the fourth pattern are the 28 participants who have both *autonomy* and *organised* values missing. We'll return to this "dfttuNA" data later.  
  
  
### 2.4 Check assumptions  
  
Linear regression has several assumptions related to the data that is to be analysed:  
  
**1. Linearity of data**: Relationship between the y (response/dependent) and x (predictor/independent) variable is linear (instead of being non-linear).  
    How to investigate: Scatterplot of y and x (clear linear association). Scatterplot of residuals vs fitted (model predicted) values of y (horizontal line without distinct pattern).  
  
**2. Independence of errors**: There is no correlation between the residuals. Residuals are the difference between the observed and model predicted values of y.  
    How to investigate: Scatterplot of residuals vs model predicted (fitted) values of y (*r* is approximately zero).  
  
**3. Normality of errors**: The residuals are normally distributed.  
    How to investigate: Normal Q-Q ("Quantile-Quantile") plot.   
  
**4. Homogeneity of errors variance (homoscedasticity)**: The variance in the residuals (of y) is constant (i.e., same) regardless of the value of predictor (x).  
    How to investigate: Scatterplot of standardized residuals vs model predicted (fitted) values of y (residuals should show no pattern across x-axis).

**5. Multicollinearity**: No strong correlations among the predictor/independent (x) variables.  
   How to investigate: Correlation matrix; Variance Inflation Factor; Durbin-Watson’s d.  
  
  
#### 2.4.1 Linearity of data  
  
First we produce diagnostic values for both predictors (autonomy, score).  
  
**Autonomy**  
  
```{r lin reg assumpions and diagnostics autonomy, message = FALSE, warning = FALSE}

# http://www.sthda.com/english/articles/39-regression-model-diagnostics/161-linear-regression-assumptions-and-diagnostics-in-r-essentials/

library(tidyverse)
library(broom)
theme_set(theme_classic())

# re-run with lm()
m1_lm <- lm(competence ~ autonomy, data = dfttu)

# see model diagnostics
m1_lm_diag <- augment(m1_lm)
head(m1_lm_diag)

```
  
head() command produces "glimpse" of the first six rows in the data:  
  
fitted: fitted (predicted) values of y (competence)  
resid: residual errors (difference between the observed and fitted values)  
hat: extreme values in the predictors (x variable here is the autonomy)  
sigma: estimated residual standard deviation when corresponding observation is dropped from model  
cooksd: Cook’s distance (influential values that may indicate an outlier)  
std.resid: outliers (or extreme values) in the outcome y variable (competence)  
  
To see the last six rows, try tail() command.  
  
  
**Score**  
  
```{r lin reg assumpions and diagnostics score, message = FALSE, warning = FALSE}

# http://www.sthda.com/english/articles/39-regression-model-diagnostics/161-linear-regression-assumptions-and-diagnostics-in-r-essentials/

library(tidyverse)
library(broom)
theme_set(theme_classic())

# re-run with lm()
m2_lm <- lm(competence ~ score, data = dfttu)

# see model diagnostics
m2_lm_diag <- augment(m2_lm)
head(m2_lm_diag)

```
  
Next we plot the observed values of the DV (competence) against the fitted regression line. First predictor is *autonomy* and then *score*.    
  
**Autonomy** as predictor  
  
```{r lin reg assumpions and diagnostics plots  autonomy 1, message = FALSE, warning = FALSE}

# http://www.sthda.com/english/articles/39-regression-model-diagnostics/161-linear-regression-assumptions-and-diagnostics-in-r-essentials/

# plot residual errors (magenta) btw obs values and the fitted reg line
ggplot(m1_lm_diag, aes(autonomy, competence)) +
  geom_point() +
  stat_smooth(method = lm, se = FALSE) +
  geom_segment(aes(xend = autonomy, yend = .fitted), color = "magenta", size = 0.3) + 
  ggtitle("Competence ~ Autonomy (original)")

```

Association between competence and autonomy looks strong and linear.  
  
**Score** as predictor  
  
```{r lin reg assumpions and diagnostics plots score 1, message = FALSE, warning = FALSE}

# http://www.sthda.com/english/articles/39-regression-model-diagnostics/161-linear-regression-assumptions-and-diagnostics-in-r-essentials/

# plot residual errors (magenta) btw obs values and the fitted reg line
ggplot(m2_lm_diag, aes(score, competence)) +
  geom_point() +
  stat_smooth(method = lm, se = FALSE) +
  geom_segment(aes(xend = score, yend = .fitted), color = "magenta", size = 0.3) + 
  ggtitle("Competence ~ score (original)")

```
  
We can see that score is not a good predictor of competence. However, there is a weak positive association.   
  
Next we look into the scatterplot of residuals vs fitted (model predicted) values of y (DV). Residual (or 'error') is the difference between the  dependent ('endogenous','response') variable's observed and fitted (estimation based on the model) values. A horizontal line without distinct patterns is an indication of a linear relationship.  
  
**Autonomy**
  
  
```{r linearity autonomy via Residuals vs Fitted plot, message = FALSE, warning = FALSE}

# run linear regression
m1_lm <- lm(competence ~ autonomy, data = dfttu)

# display first (of four) plot
plot(m1_lm,1)

```
  
We have here no distinct pattern in the residual plot.  
  
  
**Score**  
  
```{r linearity score via Residuals vs Fitted plot, message = FALSE, warning = FALSE}

# run linear regression
m2_lm <- lm(competence ~ score, data = dfttu)

# display first (of four) plot
plot(m2_lm,1)

```
  
As was with competence ~ autonomy, we have no distinct pattern in the residual plot with competence ~ score.   
  
  
#### 2.4.2 Independence of errors  
  
Plots below show that there is no correlation between the residuals (red line is horizontal and close to zero).  
  
```{r, message = FALSE, warning = FALSE}

plot(m1_lm,1)
plot(m2_lm,1)

```

  
  
  
#### 2.4.3 Normality of errors  
  
Normal Q-Q plot is widely used way to see if the residuals are normally distributed. This is the case if they (approximately) follow the straight (dashed) line.  
  
**Competence ~ autonomy**  
  
Plots below show that the competence residuals (with autonomy as a predictor) follow the line quite nicely and we assume that they are normally distributed.    
  
```{r normality of competence residuals with autonomy, message = FALSE, warning = FALSE}

plot(m1_lm,2)

```
  
Not just to rely on a plot, we calculate Shapiro-Wilk and Kolmogorov-Smirnov tests.   
  
```{r test violation of normality assumption competence and autonomy, message = FALSE, warning = FALSE}

# install.packages("olsrr")
library(olsrr)

# Shapiro-Wilk and Kolmogorov-Smirnov test values should be > .05
ols_test_normality(m1_lm) # both values are above .05

# correlations btw observed and fitted values
ols_test_correlation(m1_lm) # high correlation (r = .994) as it should be 

```
  
S-W and K-S values are above .05, so the null hypothesis that the residuals are normally distributed cannot be rejected.  
  
    
**Competence ~ score**   
  
Plots below show that the competence residuals (with score as a predictor) do not follow the line as nicely as with autonomy, so we do some further checks.  
  
```{r normality of competence residuals with score, message = FALSE, warning = FALSE}

plot(m2_lm,2)

```
  
Not just to rely on a plot, we will calculate Shapiro-Wilk and Kolmogorov-Smirnov tests.   
  
```{r test violation of normality assumption competence and score, message = FALSE, warning = FALSE}

# install.packages("olsrr")
library(olsrr)

# Shapiro-Wilk and Kolmogorov-Smirnov test value should be > .05
ols_test_normality(m2_lm) # both below .05

# correlations btw observed and fitted values
ols_test_correlation(m2_lm) # high correlation .987 as it should be 

```
  
S-W and K-S values are below .05, so the null hypothesis that the residuals are normally distributed can be rejected. However, normality assumption violation is not a big problem with larger sample sizes.  
  
  
  
#### 2.4.4 Homogeneity of errors variance    
  
**Scale-Location (or spread-Location):** Homogeneity of variance of the residuals (homoscedasticity) is met if the line is horizontal and the points are equally spread.  
- Our residuals are spread equally along the ranges of predictors (so we have no problem with heteroscedasticity). If there would be a problem, we could use a log or square root transformation of the outcome variable (y):    
  lm(log(competence) ~ autonomy, data = dfttu)  
  or  
  lm(sqrt(competence) ~ autonomy, data = dfttu)   
  
```{r homogeneity of variance autonomy, message = FALSE, warning = FALSE}

# raw data competence ~ autonomy
plot(m1_lm,3)

# square root transformed competence
m1b_lm_sqrt <- lm(sqrt(competence) ~ autonomy, data = dfttu)
# log transformed autonomy
m1c_lm_log <- lm(log(competence) ~ autonomy, data = dfttu)

plot(m1b_lm_sqrt,3)
plot(m1c_lm_log,3)


```
  
First plot shows that the residuals (of competence) are spread quite equally along the range of autonomy. Doing sqrt or log transformation to DV (competence) does not improve the situation.  
  
  
```{r homogeneity of variance score, message = FALSE, warning = FALSE}

# raw data competence ~ score
plot(m2_lm,3)

# square root transformed competence
m2b_lm_sqrt <- lm(sqrt(competence) ~ score, data = dfttu)
# log transformed autonomy
m2c_lm_log <- lm(log(competence) ~ score, data = dfttu)

plot(m2b_lm_sqrt,3)
plot(m2c_lm_log,3)

```
  
  
We also leave score untransformed.  
  
  
#### 2.4.5 Multicollinearity   
  
Next we will investigate the correlations between the variables in the model.  
  
  
```{r Pearson correlations, message = FALSE, warning = FALSE}

# create new df with only two variables
dfttu2 <- dfttu[,c("autonomy","competence","score")]

# compute correlation matrix
lavCor(dfttu2)

```
  
We see that autonomy and competence have a strong correlation (*r* = `r dec2(lavCor(dfttu2)[2,1])`).  
  
Cohen (1988) has defined some effect size "rules" for *r*:  
small     .1  
medium    .3  
large     > .5  
  
So we conclude at this point that autonomy has "a large bi-variate association" with competence. We further check if there is a multicollinearity problem with the regression models using Durbin-Watson d test.    
  
```{r, message = FALSE, warning = FALSE}

library(car)

# Null hypothesis: There is no correlation among the residuals (p-value >.05)
durbinWatsonTest(m1_lm) # com ~ aut
durbinWatsonTest(m2_lm) # com ~ sco

```
  
It seems that multicollinearity is not a problem in this data.  
  
  
#### 2.4.6 Outliers and high leverage points  
  
Final step before the linear regression analysis is to check if there are extreme values that might influence the regression results. Residuals vs Leverage plot is useful for this purpose. 
  
Outlier: extreme outcome variable (DV) value (competence).  
  - Examination of standardized residuals (residual / standard error). Potential outlier if std.res. > 3.  
  
High leverage point: extreme predictor (IV) value (autonomy, score).  
  - Examination of .hat value. Potential high leverage point if .hat > 2(p+1)/n. Note: p = number of predictors, n = number of observations.  
    - In "competence ~ autonomy" model: 2(1+1)/208 = 0.02.  
  
  
```{r outliers and high leverage points autonomy score, message = FALSE, warning = FALSE}

# plot autonomy
plot(m1_lm,5)

# plot score
plot(m2_lm,5)

```
  
First plot (competence ~ autonomy) shows that the results regarding outliers seem to be quite ok as the standardized residuals are not exceeding |3|. However, leverage values exceed 0.02 level. Top three potential outliers are participants 43, 81 and 167.   
  
Score as a predictor performs better: No outliers or high leverage point problems are present.  
  
Next we investigate the most influential values with Cook's distance method. An observation is influential if its Cook's distance value is greater than 4 / (n - p - 1). The value with our sample of 208 (n) and one predictor (p) is: 4 / (208 - 1 - 1) = 0.01941748. Also a simpler equation is applicable: 4/n = 4/208 = 0.01923077. Roughly speaking we are looking at Cook's distance values greater than 0.02.  
  
First we plot competence ~ autonomy:  
  
```{r plot Cooks distance autonomy, message = FALSE, warning = FALSE}

# autonomy

# plot 1
n <- nrow(dfttu)
plot(cooks.distance(m1_lm), main = "Cook's distance")
abline(h = 4/n, lty = 2, col = "darkred") # add cutoff line

# plot 2
plot(m1_lm, 4, 
     id.n = 10 # ask to indicate top 10 influential values
     )
abline(h = 4/n, lty = 2, col = "darkred")

```
  
Row 81 value is clearly the highest, followed by 167 and 43.  
  
And then the competence ~ score:  
  
```{r plot Cooks distance score, message = FALSE, warning = FALSE}

# score

# plot 1
n <- nrow(dfttu)
plot(cooks.distance(m2_lm), main = "Cook's distance")
abline(h = 4/n, lty = 2, col = "darkgreen")

# plot 2
plot(m2_lm, 4, 
     id.n = 10 # ask to indicate top 10 influential values
     )
abline(h = 4/n, lty = 2, col = "darkgreen")

```
  
  
Plots above show that some rows in the df have higher values than 0.02. So we proceed to further check their behavior ..    
  
  
```{r further examination of Cooks and Mahalanobis distances autonomy, message = FALSE, warning = FALSE}

# install.packages("performance")
library(performance)

# ---------------
# Cook's distance
# https://www.statology.org/how-to-identify-influential-data-points-using-cooks-distance/
# measures how much all of the fitted values in the model change when the one data point (in turn) is deleted.

# check outliers with Cook's distance, m1_lm (competence ~ autonomy)
check_outliers(m1_lm, method = "cook") # result is OK, no outliers

# display Cook's distance values in descending order
m1_lm_outlier_cook_list <- as.data.frame(check_outliers(m1_lm, method = "cook"))
m1_lm_outlier_cook_list[order(-m1_lm_outlier_cook_list$Distance_Cook),]

# we see that observation 81 has the highest value (0.058), but also several other values are above the 0.02 level
# however, the test indicates that we do not have an outlier problem here




# --------------------
# Mahalanobis distance

# to get a second opinion, we turn to Mahalanobis distance
# https://www.statology.org/mahalanobis-distance-r/
check_outliers(dfttu[c("competence","autonomy")], method = "mahalanobis")

# we have an indication of one possible outlier, participant 81
# this finding is against the Cook distance result above, but resonates with the earlier plot

# display Mahalanobis values in descending order
m1_lm_outlier_mah_list <- as.data.frame(check_outliers(dfttu[c("competence","autonomy")], method = "mahalanobis"))
m1_lm_outlier_mah_list <- m1_lm_outlier_mah_list[order(-m1_lm_outlier_mah_list$Distance_Mahalanobis),]
# add also p-values, usually we look for values smaller than .001
m1_lm_outlier_mah_list$pvalue <- pchisq(m1_lm_outlier_mah_list$Distance_Mahalanobis, 
                                        df=1, # number of variables minus one (2-1=1)
                                        lower.tail=FALSE
                                        )
m1_lm_outlier_mah_list

# p-values are not smaller than .001, so the end result is in accordance with Cook's distance: no outliers here 

```

```{r further examination of Cooks and Mahalanobis distances score, message = FALSE, warning = FALSE}

# install.packages("performance")
library(performance)

# ---------------
# Cook's distance

# check outliers with Cook's distance, m2_lm (competence ~ score)
check_outliers(m2_lm, method = "cook") # result is OK, no outliers

# display Cook distance values in descending order
m2_lm_outlier_cook_list <- as.data.frame(check_outliers(m2_lm, method = "cook"))
m2_lm_outlier_cook_list[order(-m2_lm_outlier_cook_list$Distance_Cook),]

# we see that observations 112 and 154 have the highest value (0.043), but also several other values are above the 0.02 level
# however, the test indicates that we do not have an outlier problem here


# --------------------
# Mahalanobis distance

# to the a second opinion, we turn to Mahalanobis distance
# https://www.statology.org/mahalanobis-distance-r/
check_outliers(dfttu[c("competence","score")], method = "mahalanobis")

# we have an indication of two possible outliers: 112 and 154
# this finding is against the Cook distance result above

# display Mahalanobis values in descending order
m2_lm_outlier_mah_list <- as.data.frame(check_outliers(dfttu[c("competence","score")], method = "mahalanobis"))
m2_lm_outlier_mah_list <- m2_lm_outlier_mah_list[order(-m2_lm_outlier_mah_list$Distance_Mahalanobis),]
# add also p-values, usually we look for values smaller than .001
m2_lm_outlier_mah_list$pvalue <- pchisq(m2_lm_outlier_mah_list$Distance_Mahalanobis, 
                                        df=1, 
                                        lower.tail=FALSE
                                        )
m2_lm_outlier_mah_list

# p-values are not smaller than .001, so the end result is in accordance with Cook distance: no outliers here 

```


  
## 3 Simple regression model  
  
Simple regression has only one x (or 'independent', 'predictor', 'exogenous', 'explanatory', 'manifest') variable and one y (or 'dependent', 'endogenous', 'response') variable.  
  
y = b0 + b1*x  
  
y  predicted value of the dependent/endogenous/response variable  
b0 intercept  
b1 regression coefficient  
x  value of the independent/predictor/exogenous/explanatory variable  
  
We specify a model that predicts competence with autonomy.  
  
Lectures [(PPT)](https://homepages.tuni.fi/petri.nokelainen/mmx9120/3_sem/ppt/structural_equation_modeling.ppt) contain more information about the different "roles" of the variables in SEM.  
  
  
### 3.1 Estimate the model  
  
  
```{r estimate simple regression, message = FALSE, warning = FALSE}

# specify simple regression model
regression_simple <- '
      # regression
      competence ~ 1 + autonomy
      # variance
      autonomy ~~ autonomy
      '

# estimate model with lavaan
regression_simple_fit <- sem(regression_simple,
                        data=dfttu)

```
  
```{r display summary of simple regression results, message = FALSE, warning = FALSE}

# display unstandardized results in LISREL format
summary(regression_simple_fit, 
        rsq=TRUE) # display R^2

```
  
  
In the output we see several values, of which the most important are:    
  
Estimate: Parameter value for each model parameter (here unstandardized "B")  
  
Std.err:  Standard error for each estimated parameter (average distance of observed values from the regression line)    
  
Z-value:  Wald statistic (Z = Estimate / Std.err)  
  
P(>|z|)): p-value for testing the null hypothesis that the parameter equals zero in the population  
  
R-Square: Explained variance of the dependent variable (DV)  
  
  
  
According to the output above, the intercept (average value) of autonomy is 3.786. Interpretation of the intercept of competence is a bit odd as the result now says that "Competence gets the value of 1.061 when predictor variable (autonomy) is set to zero." We do not have a zero value for autonomy (range 1-5), but we can create it with "centering grand mean" (CGM): autonomy - mean(autonomy).  
  
```{r centering autonomy, message = FALSE, warning = FALSE}

# of course we can change autonomy and competency values from 1-5 -> 0-4
dfttu$autonomy04 <- dfttu$autonomy - 1
dfttu$competence04 <- dfttu$competence - 1

# but still, it makes more sense to interpret zero as an average value of autonomy
# Note: usually we do not center the DV!

# centering by the grand mean
dfttu$autonomyCGM1 <- dfttu$autonomy - mean(dfttu$autonomy) 
# or
# install.packages("misty")
library(misty) # center()
dfttu$autonomyCGM2 <- dfttu$autonomy # we need to do this extra step as center() writes new values directly into df given in the first parameter
dfttu$autonomyCGM2 <- center(dfttu$autonomyCGM2, # function writes new values into this variable
                             type = "CGM", # center grand mean, other option is "CWC" (center withing cluster) 
                             cluster = NULL, # CWC needs the cluster name (e.g., "group")
                             value = NULL # wish to center to a specific value?
                             )

# check how different "autonomy" variants look now
# table
print(dfttu[c(1:5),c("autonomy","autonomyCGM1","autonomyCGM2")])
# plots
par(mfrow = c(1, 2)) # place plots in one row and two columns
plot(density(dfttu$autonomy), col = "black", main = "Autonomy original")
plot(density(dfttu$autonomyCGM1), col = "blue", main = "Autonomy CGM")

# original mean value of autonomy
mean(dfttu$autonomy) # 3.786058

# grand mean centered mean value of autonomy
mean(dfttu$autonomyCGM1) # -2.049622e-16 # naturally we have here exactly the same values
mean(dfttu$autonomyCGM2) # -2.049622e-16

```
  
As the plot shows, the shape of the density distribution does not change, only the average value of autonomy moves to zero.   
  
Now we can re-run the analysis ..  
    
```{r re-run simple regression, message = FALSE, warning = FALSE}

# specify simple regression
regression_simple_CGM <- '
      # regression
      competence ~ 1 + autonomyCGM1
      # variance
      autonomyCGM1 ~~ autonomyCGM1
      '

# estimate model
regression_simple_fit_CGM <- sem(regression_simple_CGM,
                        data=dfttu,
                        group = NULL, # no grouping variable (e.g., gender)
                        cluster = NULL, # no clustering variable (e.g., multilevel data may contain several values for each participant)
                        missing = "fiml.x" # full information maximum likelihood approach (good if there are missing observations)
                        )
```

  
```{r display results of the second run, message = FALSE, warning = FALSE}

# display unstandardized results in LISREL format
summary(regression_simple_fit_CGM, 
        rsq=TRUE) # display R^2

```
  
Results stay the same, but now we can say that "competence has an average value of 3.9 when autonomy is on its average level."  
  
  
```{r extract simple regression output values for inline text, message = FALSE, warning = FALSE}

# extract values from the lavaan output for the Results section

# unstandardized output values of competence ~ autonomy
reg_m1_aut_est <- parameterEstimates(regression_simple_fit_CGM)$est[2] # estimate, "B"
reg_m1_aut_se <- parameterEstimates(regression_simple_fit_CGM)$se[2] # Standard error of the estimate
reg_m1_aut_z <- parameterEstimates(regression_simple_fit_CGM)$z[2] # Z value
reg_m1_aut_p <- parameterEstimates(regression_simple_fit_CGM)$pvalue[2] # p-value
reg_m1_aut_r2 <- as.numeric(lavaan::inspect(regression_simple_fit_CGM,'r2')[1]) # R^2 value Note: need to use "as.numeric()" to extract numerical value

```
  
  
  
#### 3.4 Report the results 
  
Results of linear regression with *R lavaan* (Rosseel, 2012) showed that autonomy was statistically significant positive predictor of competence, *B* = `r dec3(reg_m1_aut_est)`, *Z* = `r dec3(reg_m1_aut_z)`, *p* = `r dec3(reg_m1_aut_p)`. This unstandardized result indicates that one unit increase in autonomy increases competence by `r dec3(reg_m1_aut_est)` units. This model explains `r dec2(reg_m1_aut_r2*100)`% of the variance of competence variable.  
  
  
#### 3.5 Missing data  
  
What if we had some missing values? We created earlier a data frame ("dfttuNA") with approximately 40% of missing values in autonomy variable. Let's re-run the analysis with that data.    
   
```{r estimate simple regression with missing data, message = FALSE, warning = FALSE}

# estimate model with lavaan
regression_simple_fit_NA <- sem(regression_simple,
                        data=dfttuNA)

```
  
```{r display summary of simple regression results with missing data, message = FALSE, warning = FALSE}

# display unstandardized results in LISREL format
summary(regression_simple_fit_NA, 
        rsq=TRUE) # display R^2

```
  
Output shows that only 129 rows (of the complete data of 208 rows) have been analysed. Results stay quite the same ... except that the standard errors are bigger (indicating worse predicting power of the model).      
  
Instead of applying listwise deletion to the rows with missing data, we can use multiple imputation (MI) for the missing values. However, it's use is not that straightforward as it generates several imputed data sets that need to be analysed at the same time (as pooled data). This requires additional steps in the modeling process.  

  
#### 3.5.1 Full information maximum likelihood  
  
Full Information Maximum Likelihood (FIML) estimator (e.g., [Cham et al., 2017](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5489914/)) is a good alternative to MI, as it is part of the model estimation process in most SEM software packages (e.g., Mplus, AMOS, lavaan).  
  
Next we utilize FIML to analyze the data. FIML estimates a likelihood function for each participant based on the observed values (so there has to be at least one observed value for each participant).   


```{r estimate simple regression again with missing data and FIML, message = FALSE, warning = FALSE}

regression_simple_fit_NA_FIML <- sem(regression_simple,
                        data=dfttuNA,
                        group = NULL, # no grouping variable (e.g., gender)
                        cluster = NULL, # no clustering variable (e.g., multilevel data may contain several values for each participant)
                        missing = "fiml", # full information maximum likelihood approach (good if there are missing observations)
                        fixed.x = FALSE # instead of using the default fixed.x = TRUE, we ask to estimate the means and variances of the variables
                        )

```

```{r display summary of simple regression results with missing data and FIML, message = FALSE, warning = FALSE}

# display unstandardized results in LISREL format
summary(regression_simple_fit_NA_FIML, 
        rsq=TRUE) # display R^2

```
  
Summary shows that there are two missing data patterns: First one is the complete observations pattern and the second one has missing values in the autonomy variable. Number of observations is again 208: FIML allows the analysis of all rows (participants), also those containing missing data. Results stay the same, but now the standard errors are smaller.   
  
  


  
## 4. Multiple regression model  
  
Multiple regression means that we add predictors to the model. We specify a model that predicts competence with autonomy and score. We use centered autonomy ("autonomyCGM") but retain score uncentered (as it has a meaningful zero value: it is plausible to get zero points from the math test).  
  
  
```{r estimate multiple regression, message = FALSE, warning = F}

# specify multiple regression
regression_multiple <- '
      # regression
      competence ~ 1 + autonomyCGM1 + score
      # variances
      autonomyCGM1 ~~ autonomyCGM1
      score ~~ score
      '

# estimate model
regression_multiple_fit <- sem(regression_multiple,
                        data=dfttu,
                        group = NULL, # no grouping variable (e.g., gender)
                        cluster = NULL, # no clustering variable (e.g., multilevel data may contain several values for each participant)
                        missing = "fiml.x" # full information maximum likelihood approach (good if there are missing observations)
                        )

```
  
Now we get a warning that the variances of autonomy and competence are quite different. We follow the instruction and look into their variances:  
  
```{r display variance information, message = FALSE, warning = FALSE}

# check variances
varTable(regression_multiple_fit)

# extract variances of autonomy and score
reg_m2_aut_s2 <- varTable(regression_multiple_fit)$var[2] # s2 is abbreviation of s^2 = sigma^2 = variance
reg_m2_score_s2 <- varTable(regression_multiple_fit)$var[3]

```
  
Indeed, the "var" column of the output indicates that the variance of autonomy is `r dec3(reg_m2_aut_s2)` and the variance of score is `r dec3(reg_m2_score_s2)` - quite a difference! Next we will use *R* base library's *scale()* function to center and scale both autonomy and score (i.e., set their mean *M* to zero and standard deviation *SD* to one).  
  
  
```{r autonomy score density, message = FALSE, warning = FALSE}

# We start by plotting density of autonomy and score before any centering or scaling actions

# The probability density function describes the probability of the variable taking certain value
# It is a smoothed (Parzen–Rosenblatt / kernel estimator) version of the histogram
# Default kernel is "Gaussian", but there are several options
# Bandwidth controls the smoothing: Too big values oversmooth and too small values undersmooth (overfit)
# density() uses "rule of thumb" as default method
# https://r-coder.com/density-plot-r/

# AUTONOMY
# create a histogram
hist(dfttu$autonomy, freq = F, main = "Autonomy histogram and density")
# add density line to the histogram
lines(density(dfttu$autonomy), lwd = 2, col = "blue")

# compare bandwidth settings
par(mfrow = c(1, 3)) # place plots in one row and three columns
plot(density(dfttu$autonomy, 
             kernel="gaussian", # type "?density()" in Console to see other options
             bw="nrd0"), # default bw is Rule of Thumb (in this case 0.1155)
     col = "black", main = "Aut bw=RoT")
plot(density(dfttu$autonomy, kernel="gaussian", bw=0.05),  # bw set too small
     col = "blue", main = "Aut bw=0.05")
plot(density(dfttu$autonomy, kernel="gaussian", bw=0.5), # set too big
     col = "red", main = "Aut bw=0.5")
# # other bw options if you wish to try
# plot(density(dfttu$autonomy, kernel="gaussian", bw="UCV"), # Unbiased cross validation
#      col = "green", main = "Aut bw=UCV")
# plot(density(dfttu$autonomy, kernel="gaussian", bw="SJ"), # Plug-in by Sheather and Jones
#      col = "orange", main = "Aut bw=SJ")

# using scale() function ... AND DOING NOTHING :-D
aut_ori <- scale(dfttu$autonomy, 
                center = F, # DO NOT center (F = FALSE)
                scale = F   # DO NOT scale (F = FALSE)
                )

# -----
# SCORE
# histogram + density line
hist(dfttu$score, freq = F, main = "Score histogram and density")
lines(density(dfttu$score), lwd = 2, col = "darkred")

# AUTONOMY AND SCORE
d_a <- density(dfttu$autonomy)
d_s <- density(dfttu$score)

plot(d_a, lwd = 2, col = "darkblue",
    main = "Autonomy and score densities", xlab = "",
    xlim = c(min(d_a$x, d_s$x), c(max(d_a$x, d_s$x))),  # Min and Max X-axis limits
    ylim = c(min(d_a$y, d_s$y), c(max(d_a$y, d_s$y))))  # Min and Max Y-axis limits

lines(d_s, col = "darkred", lwd = 2)

```
  
Last figure clearly shows that there is a need for centering and scaling if we use autonomy and score in the same analysis as predictors.  
  
  
```{r center autonomy, message = FALSE, warning = FALSE}

# -----------------
# CENTERED autonomy
# center = subtract mean of autonomy from each autonomy value
# to produce M = 0

# equation approach
dfttu$autonomyCentered <- dfttu$autonomy - mean(dfttu$autonomy)

# using scale() function
aut_cen <- scale(dfttu$autonomy, 
                center = T, # center (T = TRUE)
                scale = F   # do not scale (F = FALSE)
                )

# plot centered autonomy
plot(density(aut_cen), col = "black", main = "Autonomy Centered")
```
  
Centering works as it should: The shape of distribution is the same, but the average value has moved to zero.  
  

```{r scale autonomy, message = FALSE, warning = FALSE}

# ---------------
# SCALED autonomy
# scale = divide each autonomy value by 
# root mean square RMS (scale function default if not centered)
# or
# standard deviation SD (scale function default if centered)
# to make all variables comparable (by one unit change)

# equation approach
# RMS
dfttu$autonomyScaledRMS <- dfttu$autonomy / sqrt(sum(dfttu$autonomy^2)/(length(dfttu$autonomy)-1))
# SD
dfttu$autonomyScaledSD <- dfttu$autonomy / sd(dfttu$autonomy)

# using scale() function
# RMS
aut_scaRMS <- scale(dfttu$autonomy, 
                 center = F, # do not center (F = FALSE)
                 scale = T   # scale (T = TRUE)
                 )
# SD
aut_scaSD <- scale(dfttu$autonomy, 
                 center = F, 
                 scale = apply(dfttu[,4,drop=F], # fourth var (autonomy) in df, dimensions of df are kept
                               2,  # 1=apply to rows, 2=apply to columns
                               sd) # calculate SD 
                 )

# plot RMS and SD scaled autonomy
d_a_scaRMS <- density(aut_scaRMS)
d_a_scaSD <- density(aut_scaSD)

plot(d_a_scaRMS, lwd = 2, col = "darkblue",
    main = "Autonomy RMS (blue) and SD (red) scaled densities", xlab = "",
    xlim = c(min(d_a_scaRMS$x, d_a_scaSD$x), c(max(d_a_scaRMS$x, d_a_scaSD$x))),  # Min and Max X-axis limits
    ylim = c(min(d_a_scaRMS$y, d_a_scaSD$y), c(max(d_a_scaRMS$y, d_a_scaSD$y))))  # Min and Max Y-axis limits

lines(d_a_scaSD, col = "darkred", lwd = 2)

```


```{r center and scale autonomy, message = FALSE, warning = FALSE}

# ----------------------------
# CENTERED and SCALED autonomy

# equation approach
dfttu$autonomyCenSca <- (dfttu$autonomy - mean(dfttu$autonomy)) / sd(dfttu$autonomy)

# using scale() function
aut_cen_sca <- scale(dfttu$autonomy, center = T, scale = T)

# plot centered and scaled autonomy
plot(density(aut_cen_sca), col = "blue", main = "Aut Scale RMS")
```

```{r compare plots and statistics autonomy, message = FALSE, warning = FALSE}

# --------------------------------------------------------------------------
# plot 1) original, 2) centered, 3) scaled and 4) centered & scaled autonomy
par(mfrow = c(2, 2))
plot(density(dfttu$autonomy), col = "black", main = "Autonomy")
plot(density(dfttu$autonomyCentered), col = "blue", main = "Autonomy Centered")
plot(density(dfttu$autonomyScaledRMS), col = "red", main = "Autonomy Scaled")
plot(density(dfttu$autonomyCenSca), col = "green", main = "Autonomy Centered and Scaled")

# --------------------
# compare descriptives
df_scale <- data.frame(aut_ori, aut_cen, aut_scaRMS, aut_scaSD, aut_cen_sca)
library(psych)
describe(df_scale)

```

  
Figures above show that the distributional shape is the same, but the average point varies. The last (green) density line represents both centered and scaled autonomy with *M* = 0 and *SD* = 1 (this can also be seen in the last row of the table "aut_cen_sca"). Note that the *[mad](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/mad)* stands for "the median of the absolute deviations from the median".     
  
Next we do centering and scaling to score.  

```{r center and scale score, message = FALSE, warning = FALSE}

# center and scale
dfttu$scoreCenSca <- (dfttu$score - mean(dfttu$score)) / sd(dfttu$score)

# plot against the original score
d_s <- density(dfttu$score)
d_s_cenSca <- density(dfttu$scoreCenSca)

plot(d_s, lwd = 2, col = "black",
    main = "Score original (black), centered and scaled (red) densities", xlab = "",
    xlim = c(min(d_s$x, d_s_cenSca$x), c(max(d_s$x, d_s_cenSca$x))),  # Min and Max X-axis limits
    ylim = c(min(d_s$y, d_s_cenSca$y), c(max(d_s$y, d_s_cenSca$y))))  # Min and Max Y-axis limits

lines(d_s_cenSca, col = "darkred", lwd = 2)

```
  
Red density line is the centered and scaled mathematics score. Next we plot autonomy and score into same linegraph.  
  
  
```{r plot centered and scaled autonomy and score, message = FALSE, warning = FALSE}

# we run everything once more ...

# center and scale
dfttu$autonomyCenSca <- (dfttu$autonomy - mean(dfttu$autonomy)) / sd(dfttu$autonomy)
dfttu$scoreCenSca <- (dfttu$score - mean(dfttu$score)) / sd(dfttu$score)

# plot densities
da <- density(dfttu$autonomyCenSca)
ds <- density(dfttu$scoreCenSca)

plot(da, lwd = 2, col = "blue",
    main = "Centered and scaled densities of Autonomy (blue) and Score (red)", xlab = "",
    xlim = c(min(da$x, ds$x), c(max(da$x, ds$x))),  # Min and Max X-axis limits
    ylim = c(min(da$y, ds$y), c(max(da$y, ds$y)))  # Min and Max Y-axis limits
    )

lines(ds, col = "darkred", lwd = 2)

```
  
Now we have both predictors centered and scaled (mean is zero and standard deviation is one). We re-run the analysis with the new variables.     
  
  
```{r rerun analysis with scaled predictors, message = F, warning = F}

# specify multiple regression2
regression_multiple2 <- '
      # regression
      competence ~ 1 + autonomyCenSca + scoreCenSca
      # variances
      autonomyCenSca ~~ autonomyCenSca
      scoreCenSca ~~ scoreCenSca
      '

# estimate model2
regression_multiple_fit2 <- sem(regression_multiple2,
                        data=dfttu,
                        group = NULL, # no grouping variable (e.g., gender)
                        cluster = NULL, # no clustering variable (e.g., multilevel data may contain several values for each participant)
                        missing = "fiml.x" # full information maximum likelihood approach (good if there are missing observations)
                        )

```

```{r check again variances, message = FALSE, warning = FALSE}

# check variances
varTable(regression_multiple_fit2)

```
  
Now the mean for both scaled predictors is zero and the sd is one. Let's see the outputs of the original (uncentered/scaled) and new analyses ..    
  
```{r display summary multiple regression, message = FALSE, warning = FALSE}

# display unstandardized results in LISREL format
# original variables
summary(regression_multiple_fit, 
        rsq=TRUE) # display R^2

# scaled variables
summary(regression_multiple_fit2, 
        rsq=TRUE) # display R^2

```
  
Results are identical and the interpretation stays the same - except that after centering and scaling the intercept value of competence refers to a situation where predictors (autonomy, score) are on their average level (*M* = 0).  
  
```{r  display unstandardised results multiple regression, message = FALSE, warning = FALSE}

# unstandardized ("B")
# uncentered and unscaled (= original)
parameterEstimates(regression_multiple_fit)
# centered and scaled
parameterEstimates(regression_multiple_fit2)

```
  
Although unstandardized estimates differ after centering and scaling, conclusions remain similar: Autonomy is a stat sig predictor of competence.  
  
  
```{r  display standardised results multiple regression, message = FALSE, warning = FALSE}

# standardized ("beta")
# uncentered and unscaled (= original)
standardizedSolution(regression_multiple_fit)
# centered and scaled
standardizedSolution(regression_multiple_fit2)

```

  
Standardized estimates between original and centered/scaled predictors are (not surprisingly) quite the same, as the standardization is also applied to the results of original (uncentered and unscaled) variables.  
    

```{r store unstandardized and standardized values multiple regression, message = FALSE, warning = FALSE}

# store predictor estimates for inline text
# https://lavaan.ugent.be/tutorial/inspect.html

#----------------------
# unstandardized values

# autonomy 
# estimate
reg_m2_aut_est <-  parameterEstimates(regression_multiple_fit)$est[2] # stores estimate 0.7530086
# or 
reg_m2_aut_est <- regression_multiple_fit@ParTable$est[2]

# standard error
reg_m2_aut_se <-  parameterEstimates(regression_multiple_fit)$se[2] # stores standard error 0.045
# or
reg_m2_aut_se <- regression_multiple_fit@ParTable$se[2]

# p-value
reg_m2_aut_p <-  parameterEstimates(regression_multiple_fit)$pvalue[2] # stores p-value < . 001


# score 
# estimate
reg_m2_score_est <-  parameterEstimates(regression_multiple_fit)$est[3] # 0.001587065

# standard error
reg_m2_score_se <-  parameterEstimates(regression_multiple_fit)$se[3] # 0.0009398335

# p-value
reg_m2_score_p <-  parameterEstimates(regression_multiple_fit)$pvalue[3] # 0.0912835


#--------------------
# standardized values
# these are usually presented in journal articles

# autonomy 
# estimate
reg_m2_aut_est.std <-  standardizedSolution(regression_multiple_fit)$est.std[2] # stores estimate 0.7546165

# standard error
reg_m2_aut_se.std <-  standardizedSolution(regression_multiple_fit)$se[2] # stores standard error 0.02974363

# p-value
reg_m2_aut_p.std <-  standardizedSolution(regression_multiple_fit)$pvalue[2] # stores p-value < . 001

# score estimate
reg_m2_score_est.std <-  standardizedSolution(regression_multiple_fit)$est.std[3] # 0.07631541

# score standard error
reg_m2_score_se.std <-  standardizedSolution(regression_multiple_fit)$se[3] # 0.04519976

# score p-value
reg_m2_score_p.std <-  standardizedSolution(regression_multiple_fit)$pvalue[3] # 0.09133388


#_-----------------------
# R squared for the model
# percentage of the dependent variable variation that the (linear) model explains
reg_m2_r2 <- as.numeric(lavaan::inspect(regression_multiple_fit,'r2')[1])

```
  
### 4.1 Report the results     
  
Results of multiple regression analysis (Figure 1) with *R lavaan* showed that autonomy was a significant predictor of competence, $\beta$ = `r dec3(reg_m2_aut_est.std)`, *p* = `r dec3(reg_m2_aut_p.std)`. According to the model, one standard deviation increase in autonomy (*SD* = `r dec2(sd(dfttu$autonomy))`) produces `r dec3(reg_m2_aut_est.std)` SD increase in competence (*SD* = `r dec2(sd(dfttu$competence))`). In practice this means that when a student's autonomy level goes up by `r dec2(sd(dfttu$autonomy))` units (on a 5-point self-assessment scale), then his/her competence goes up by `r dec3(reg_m2_aut_est.std)` * `r dec2(sd(dfttu$competence))` = `r dec3(reg_m2_aut_est.std * sd(dfttu$competence))` units (on a similar 5-point self-assessment scale). Although the mathematics score was a non significant positive predictor of competence, we can make a similar prediction on its effect: When a student's mathematics score goes up by 1 *SD* (`r dec2(sd(dfttu$score))` points (on a scale from 0-100), his/her competence goes up by `r dec3(reg_m2_score_est.std)` * `r dec2(sd(dfttu$competence))` = `r dec3(reg_m2_score_est.std * sd(dfttu$competence))` units (on a 5-point self-assessment scale).   
  
Standard errors (average distance of observed values from the regression line) were quite small for autonomy (`r dec3(reg_m2_aut_se.std)`) and mathematics score (`r dec3(reg_m2_score_se.std)`). Autonomy and mathematics score explained together `r dec2(reg_m2_r2*100)`% of the variance of competence. The improvement was quite small compared to the previous model that had only one predictor (autonomy, *R<sup>2</sup>* = `r dec2(reg_m1_aut_r2*100)`%).  
  
  
Figure 1. Multiple regression of competence predicted by autonomy and mathematics score (standardized estimates).  
   
```{r diagram multiple regression, message = FALSE, warning = FALSE}

# produce path diagram
semPaths(regression_multiple_fit,
         layout="tree",
         "std", # unstandardized parameter estimates OR use "par" for standardized parameter estimates
         # "par", 
         edge.label.cex=0.7,
         edge.width=0.2,
         edge.color="black",
         # optimizeLatRes=TRUE,
         label.scale=FALSE, # do not make labels smaller to fit inside boxes
         # normalize=TRUE,
         fade=FALSE, # do not fade non sig lines
         as.expression=c("nodes"), # display full labels instead of abbreviated
         # levels=c(1,3,6,9) # set distances between elements of the graph
         )

```
  
## Exercise 1  
  
Report here the **unstandardized** results of the above analysis: Describe the "B" estimates and produce a diagram/figure of the model with unstandardized values.    
  
  
## 5. Multivariate regression model  
  
Multivariate regression model has more than one endogenous variable (DV). We specify a model that predicts both competence (endogenous, DV) and autonomy (endogenous, DV)  with deep approach to learning (exogenous, IV) and organised approach to studying (exogenous, IV). Now that autonomy is not a predictor anymore, we use a non-centered version of it (*autonomy*).    
  
  
  
```{r estimate multivariate regression, message = FALSE, warning = FALSE}

# center deep and organised by grand mean
dfttu$deepCGM <- dfttu$deep - mean(dfttu$deep) 
dfttu$organisedCGM <- dfttu$organised - mean(dfttu$organised) 

# specify multivariate regression
regression_multivariate <- '
 # regressions
   autonomy ~ 1 + deepCGM
   competence ~ 1 + deepCGM + organisedCGM
   
 # variances and covariance
   deepCGM ~~ deepCGM
   organisedCGM ~~ organisedCGM
   deepCGM ~~ organisedCGM
'

# estimate model
regression_multivariate_fit <- sem(regression_multivariate, 
                                   data=dfttu,
                                   group = NULL, # no grouping variable (e.g., gender)
                                   cluster = NULL, # no clustering variable (e.g., multilevel data may contain several values for each participant)
                                   missing = "fiml.x" # full information maximum likelihood approach (good if there are missing observations)
                        )

```
  
```{r display summary multivariate regression, message = FALSE, warning = FALSE}

# display results in LISREL format
summary(regression_multivariate_fit, 
        rsq=TRUE) # display R^2

```

```{r  display un/standardized results multivariate regression, message = FALSE, warning = FALSE}

# unstandardized
parameterEstimates(regression_multivariate_fit)

# standardized
standardizedSolution(regression_multivariate_fit)

```


```{r standardized values multivariate regression, message = FALSE, warning = FALSE}

#--------------------
# standardized values

# aut ~ deep estimate
reg_m3_aut_deep_est.std <-  standardizedSolution(regression_multivariate_fit)$est.std[2]

# aut ~ deep standard error
reg_m3_aut_deep_se.std <-  standardizedSolution(regression_multivariate_fit)$se[2]

# aut ~ deep p-value
reg_m3_aut_deep_p.std <-  standardizedSolution(regression_multivariate_fit)$pvalue[2]


# com ~ deep estimate
reg_m3_com_deep_est.std <-  standardizedSolution(regression_multivariate_fit)$est.std[4]

# com ~ deep standard error
reg_m3_com_deep_se.std <-  standardizedSolution(regression_multivariate_fit)$se[4]

# com ~ deep p-value
reg_m3_com_deep_p.std <-  standardizedSolution(regression_multivariate_fit)$pvalue[4]


#_-----------------------
# R squared for the model
# percentage of the dependent variable variation that the (linear) model explains
reg_m3_aut_r2 <- lavaan::inspect(regression_multivariate_fit,'r2')[1]
reg_m3_com_r2 <- lavaan::inspect(regression_multivariate_fit,'r2')[2]


```
  
### 5.1 Report the results     
  
Results of multivariate regression analysis (Figure 2) with *R lavaan* showed that deep approach to learning was a significant positive predictor for both autonomy ($\beta$ = `r dec3(reg_m3_aut_deep_est.std)`, *p* = `r dec3(reg_m3_aut_deep_p.std)`) and competence ($\beta$ = `r dec3(reg_m3_com_deep_est.std)`, *p* = `r dec3(reg_m3_com_deep_p.std)`). Organised approach to studying ... **[Add results here]**. Deep approach to studying explained **[Add]**% of the variance of autonomy. Deep and organised explained together **[Add]**% of the variance of competence.    
  
  
## Exercise 2    
  
Add "missing" results to the above analysis.  
  
  
Figure 2. Multivariate regression of competence and autonomy predicted by deep approach to learning and organised approach to studying (standardized estimates).  
   
```{r diagram multivariate regression, message = FALSE, warning = FALSE}

# produce diagram
semPaths(regression_multivariate_fit,
         layout="tree",
         "std", # unstandardized parameter estimates OR "par" for standardized parameter estimates
         edge.label.cex=0.7,
         edge.width=0.2,
         edge.color="black",
         # optimizeLatRes=TRUE,
         label.scale=FALSE, # do not make labels smaller to fit inside boxes
         # normalize=TRUE,
         fade=FALSE, # do not fade non sig lines
         as.expression=c("nodes"), # display full labels instead of abbreviated
         # levels=c(1,3,6,9) # set distances between elements of the graph
         )

```
  
  
## 6. Path model  
  
Multivariate regression model becomes a path model after endogenous variables are allowed to predict other endogenous variables. We specify a model that predicts 1) competence with deep approach to learning, organised approach to studying and autonomy (which is the "path" part of the model), and 2) autonomy with deep approach to learning.  
  
```{r produce variance-covariance matrix, message = FALSE, warning = FALSE}

# create new df with the four variables in the path model
dfttu3 <- dfttu[,c("autonomy","competence","deep", "organised")]

# produce variance-covariance matrix
cov(dfttu3)

# store some interesting values for the inline text

# mean of autonomy
dfttu_m_aut <- mean(dfttu3$autonomy)

# variance of autonomy
dfttu_var_aut <- cov(dfttu3)[1,1]

# square root of the variance of autonomy (SD)
dfttu_sd_aut <- sqrt(cov(dfttu3)[1,1])

```
  
Diagonal values represent variances and other values covariances. Variance-covariance matrix is symmetric (same values below and above the diagonal). Square root of the diagonal values produce standard deviations. For example, autonomy variance ($\sigma^2$) is `r dec3(dfttu_var_aut)` and *SD* ($\sigma$) is `r dec2(dfttu_sd_aut)`. As the mean of the autonomy is `r dec2(dfttu_m_aut)`, +/- 1*SD* gives range where approximately 68% of the autonomy values are. All the covariances are positive (e.g., if autonomy increases, also competence increases).  
  
  
```{r estimate path model, message = FALSE, warning = FALSE}

# specify regressions
path_model <- '
   autonomy ~ 1 + deepCGM
   competence ~ 1 + deepCGM + organisedCGM + autonomy # we add autonomy here as a predictor of competence
   
 # variances and covariances
   deepCGM ~~ deepCGM
   organisedCGM ~~ organisedCGM
   deepCGM ~~ organisedCGM
'

# estimate model
path_model_fit <- sem(path_model, 
                      data = dfttu,
                      group = NULL, # no grouping variable (e.g., gender, control/intervention,...)
                      cluster = NULL, # no clustering variable (e.g., multilevel data may contain several values for each participant)
                      missing = "fiml.x" # full information maximum likelihood approach (good if there are missing observations)
                      )

```
  
```{r display summary path model, message = FALSE, warning = FALSE}

# display results in LISREL format
summary(path_model_fit, 
        rsq=TRUE) # display R^2

```

```{r  display un/standardized results path model, message = FALSE, warning = FALSE}

# unstandardized
parameterEstimates(path_model_fit)

# standardized
standardizedSolution(path_model_fit)

```
  
```{r unstandardized values path model, message = FALSE, warning = FALSE}

#--------------------
# unstandardized values

# first, we investigate the structure of std output object
str(parameterEstimates(path_model_fit)) # we learn that "$est.sdt", "$se" and "$pvalue" contain estimates that we are interested in reporting

# second, we look into the std output to see the "rows" we wish to refer
parameterEstimates(path_model_fit) 

# third, we store values from the output into variables for inline text referencing ...

# autonomy is only predicted by deep
# aut ~ deep 
# estimate
path_m4_aut_deep_est <-  parameterEstimates(path_model_fit)$est[2] # row 2 in the column "est.sdt" contains beta (standardized) estimate of autonomy ~ deep approach to learning
# standard error
path_m4_aut_deep_se <-  parameterEstimates(path_model_fit)$se[2]
# p-value
path_m4_aut_deep_p <-  parameterEstimates(path_model_fit)$pvalue[2]


# competence is predicted by deep, organised and autonomy
# com ~ deep
# estimate
path_m4_com_deep_est <-  parameterEstimates(path_model_fit)$est[4]
# standard error
path_m4_com_deep_se <-  parameterEstimates(path_model_fit)$se[4]
# p-value
path_m4_com_deep_p <-  parameterEstimates(path_model_fit)$pvalue[4]

# com ~ organised
# estimate
path_m4_com_org_est <-  parameterEstimates(path_model_fit)$est[5]
# standard error
path_m4_com_org_se <-  parameterEstimates(path_model_fit)$se[5]
# p-value
path_m4_com_org_p <-  parameterEstimates(path_model_fit)$pvalue[5]

# com ~ autonomy
# estimate
path_m4_com_aut_est <-  parameterEstimates(path_model_fit)$est[6]
# standard error
path_m4_com_aut_se <-  parameterEstimates(path_model_fit)$se[6]
# p-value
path_m4_com_aut_p <-  parameterEstimates(path_model_fit)$pvalue[6]


# -----------------------
# R squared for the model
# percentage of the dependent variable variation that the (linear) model explains
path_m4_aut_r2 <- lavaan::inspect(path_model_fit,'r2')[1] # how much "deep" explains the variance of "autonomy"?
path_m4_com_r2 <- lavaan::inspect(path_model_fit,'r2')[2] # how much "deep", "organised" and "autonomy" explain the variance of "competence"?

```
  
```{r standardized values path model, message = FALSE, warning = FALSE}

#--------------------
# standardized values

# first, we investigate the structure of std output object
str(standardizedSolution(path_model_fit)) # "$est.sdt", "$se" and "$pvalue" contain estimates that we are interested in reporting

# second, we look into the std output to see the "rows" we wish to refer
standardizedSolution(path_model_fit) 

# third, we store values from the std output into variables for inline text referencing

# aut ~ deep 
# estimate (row 2 in the column "est.sdt" contains beta (standardized) estimate of autonomy ~ deep approach to learning)
path_m4_aut_deep_est.std <-  standardizedSolution(path_model_fit)$est.std[2] 
# standard error
path_m4_aut_deep_se.std <-  standardizedSolution(path_model_fit)$se[2]
# p-value
path_m4_aut_deep_p.std <-  standardizedSolution(path_model_fit)$pvalue[2]

# com ~ deep 
# estimate
path_m4_com_deep_score_est.std <-  standardizedSolution(path_model_fit)$est.std[4]
# standard error
path_m4_com_deep_se.std <-  standardizedSolution(path_model_fit)$se[4]
# p-value
path_m4_com_deep_p.std <-  standardizedSolution(path_model_fit)$pvalue[4]

# com ~ organised
# estimate
path_m4_com_org_est.std <-  standardizedSolution(path_model_fit)$est.std[5]
# standard error
path_m4_com_org_se.std <-  standardizedSolution(path_model_fit)$se[5]
# p-value
path_m4_com_org_p.std <-  standardizedSolution(path_model_fit)$pvalue[5]

# com ~ autonomy
# estimate
path_m4_com_aut_est.std <-  standardizedSolution(path_model_fit)$est.std[6]
# standard error
path_m4_com_aut_se.std <-  standardizedSolution(path_model_fit)$se[6]
# p-value
path_m4_com_aut_p.std <-  standardizedSolution(path_model_fit)$pvalue[6]

```
  
### 6.1 Report the results    
  
  
Results of path analysis (Figure 3) with *R lavaan* showed that deep approach to learning was a positive statistically significant predictor of autonomy,  *B* = `r dec3(path_m4_aut_deep_est)`, *SE* = `r dec3(path_m4_aut_deep_se)`, *p* = `r dec3(path_m4_aut_deep_p)`. This unstandardized result indicates that a one unit increase in deep increases autonomy by `r dec3(path_m4_aut_deep_est)` units. Deep approach to learning explained `r dec1(path_m4_aut_r2*100)`% of the variance of autonomy variable. It was also investigated if deep approach to learning, organised approach to studying and autonomy are related to competence. Results showed that all three predictors had positive associations with competence, but only deep (*B* = `r dec3(path_m4_com_deep_est)`, *SE* = `r dec3(path_m4_com_deep_se)`, *p* = `r dec3(path_m4_com_deep_p)`) and autonomy (*B* = `r dec3(path_m4_com_aut_est)`, *SE* = `r dec3(path_m4_com_aut_se)`, *p* = `r dec3(path_m4_com_aut_p)`) were statistically significant. These variables explained `r dec2(path_m4_com_r2*100)`% of the variance of competence variable.  
  
    
Figure 3. Path model of competence and autonomy with selected predictors (unstandardized estimates).  
   
```{r diagram path model, message = FALSE, warning = FALSE}

# produce diagram
semPaths(path_model_fit,
         layout="tree",
         "par", # unstandardized parameter estimates OR "par" for standardized parameter estimates
         edge.label.cex=0.7,
         edge.width=0.2,
         edge.color="black",
         # optimizeLatRes=TRUE,
         label.scale=FALSE, # do not make labels smaller to fit inside boxes
         # normalize=TRUE,
         fade=FALSE, # do not fade non sig lines
         as.expression=c("nodes"), # display full labels instead of abbreviated
         # levels=c(1,3,6,9) # set distances between elements of the graph
         )

```

  
## Exercise 3  
  
Produce text that describes the standardized results of the analyses above (alongside with a figure).  
  
  
  
## 7. SEM model  
  

```{r read new data, message = FALSE, warning = FALSE}

dfttu4 <- read.csv("data1_csv.dat", header = T)
summary(dfttu4)

```
  
Our data is from a Finnish polytechnics of higher education, collected in 2000 (n = 447). The respondents are the staff members of the organization. The original measurement instrument (Growth-oriented Atmosphere Questionnaire, GOAQ) has 13 factors and 92-items (Ruohotie, 1996; Ruohotie, Nokelainen & Tirri, 2002), but for the purposes of this exercise we have selected the following factors and sample items. Five –point scale from 1 (totally disagree) to 5 (totally agree) was used in the online survey.  
  
  
**Factor 1: Psychic stress of the work (PSY)**     
v44 I feel that I am beginning to dislike my work.     
v45 I feel that it is getting more difficult for me to take the initiative.     
  
**Factor 2: Encouraging leadership (ENC)**     
v5  My manager is friendly and easily approachable.     
v7  My manager pays attention to my suggestions and wishes.     
v17 My manager works with a team to find solutions.     
  
**Factor 3: Build-up of work requirements (BUI)**     
v42r My workload has increased during the past years.     
v43  My working pace has increased in recent years.   
  
  
```{r correlations sem model, message = FALSE, warning = FALSE}

# create new df with only relevant variables for correlation matrix
dfttu4cor <- dfttu4[,c("v44","v45","v5","v7","v17","v42r","v43")]

# compute correlation matrix
lavCor(dfttu4cor)

```

Correlation matrix shows that the observed variables for each factor have quite strong and positive correlations, as expected. For example, "v44" and "v45" are manifestations of the same latent variable ("PSY") and so they have a strong positive correlation (*r* = `r dec3(lavCor(dfttu4cor)[1,2])`). Also the negative correlations make sense: It is plausible that the three observed variables ("v5","v7","v17") that belong to the latent ENC variable (Encouraging leadership) are negatively associated to observed variables related to psychic stress ("v44","v45") and build-up of work requirements ("v42r","v43").  
  
However, there is one suspicious correlation that we need to examine further: "v44" (PSY) and "v43" (BUI) are quite strongly positively correlated ( *r* = `r dec3(lavCor(dfttu4cor)[1,7])`) although they belong to different factors?! We will return to this issue in the last section of this exercise (via "modification indices").  
  
 
### 7.1 Specify the SEM model  
  
  
According to Beaujean (2014, p. 37), "A structural equation model (SEM) is a broad class of statistical models that consists of two parts: the structural model and the latent variable model. The structural model consists of the regression-like relationships among the variables (i.e., the path analysis). The latent variable model (LVM) forms the latent variables (LVs) used in the structural model. When a LVM is analyzed without a structural model, it is sometimes called a confirmatory factor analysis (CFA). If there was not a hypothesized structure for the latent variable model, then it would be an exploratory factor analysis (EFA)."  
  
The SEM model (Figure 4) aims to analyse if Encouraging leadership (ENC) and Build-up of work requirements (BUI) are related to the Psychic stress of the work (PSY).  
  
  
Figure 4. SEM model between Encouraging leadership, Build-up of work requirements and Psychic stress of the work
  
```{r sem model, message = F, echo = F}

sem_model <- '
ENC =~ NA*v5 + v7 + 1*v17   # Encouraging leadership
BUI =~ 1*v42r + v43         # Build-up of work requirements
PSY =~ 1*v44 + v45          # Psychic stress of the work
PSY ~ ENC + BUI             # Predict PSY with ENC and BUI
ENC ~~ BUI                  # Let ENC and BUI correlate
'

sem_model_fit <- sem(sem_model, 
                    data=dfttu4,
                    group = NULL, # no grouping variable (e.g., gender, control/intervention,...)
                    cluster = NULL, # no clustering variable (e.g., multilevel data may contain several values for each participant)
                    missing = "fiml.x" # full information maximum likelihood approach
                    )

semPaths(sem_model_fit, 
         layout="tree2", #spring
         as.expression=c("nodes","edges"),
         #levels=c(1,33,66,99),
         #title=TRUE,
         #title.color="black",
         sizeMan = 7, 
         sizeInt = 5,
         sizeLat = 16, 
         label.prop=0.7, 
         curve=1.5, 
         bg="white",
         groups="latents",
         pastel=TRUE,
         intercepts=FALSE,
         label.cex=1,
         label.color="black",
         label.scale=TRUE,
         borders=TRUE,
         #font=5,
         arrows=TRUE,
         optimizeLatRes=TRUE,
         #label.norm="1",
         #reorder=TRUE,
         #exoVar=TRUE,
         edge.label.cex=0.1,
         edge.width=1.5,
         width=7,
         height=7,
         normalize=TRUE
         )

```
  
  
ENC is a latent exogenous (IV) variable that is constructed with three observed variables (v5, v7, v17). The code below shows that the highest loading variable is not the first one (v5), but the last one (v17) instead. To make the model identifiable, we fix one indicator path to 1 [[see more details](http://www-personal.umich.edu/~gonzo/papers/sem.pdf)] for each latent variable. In the case of PSY, we fix the strongest loading variable (v17) instead of the first variable in the list (using "NA*" removes the automatic fix to 1 for v5).  
  
BUI is also an exogenous latent variable, but based only on two observed variables (v42r, v43). The "r" in v42r stands for a reverse coded variable: The original statement in v42 was "My workload has decreased during the past years.", but as we wish to measure build-up of work, we use reversed wording "My workload has increased during the past years." instead (and indicate that by naming the variable "v42r").  
  
PSY is an endogenous (DV) latent variable with two observed indicators (v44, v45).  
  
  
```{r specify the sem model, message = FALSE, warning = FALSE}

sem_model <- '
# Factors
ENC =~ NA*v5 + v7 + 1*v17   # Encouraging leadership
BUI =~ 1*v42r + v43         # Build-up of work requirements
PSY =~ 1*v44 + v45          # Psychic stress of the work
# Regressions
PSY ~ ENC + BUI             # Predict PSY with ENC and BUI
# Covariances
ENC ~~ BUI                  # Let ENC and BUI correlate
'

```

  
### 7.2 Estimate the model  
  
The model specified in the previous section is now estimated. The default estimation method is Maximum Likelihood (ML). As there are some missing data, we use [Full Information Maximum Likelihood (FIML)](https://stat.utexas.edu/software-faqs/general#missingincompletedata) instead of casewise deletion.  
  

```{r estimate the sem model, message = FALSE, warning = FALSE}

sem_model_fit <- sem(sem_model, 
                    data=dfttu4,
                    
                    group = NULL, # no grouping variable (e.g., gender, control/intervention,...)
                    cluster = NULL, # no clustering variable (e.g., multilevel data may contain several values for each participant)
                    missing = "fiml.x" # full information maximum likelihood approach
                    )

```
  
```{r sneak peak into results, message = FALSE, warning = FALSE}

# summary(sem_model_fit, 
#         standardized=TRUE)
# 
# # Latent Variables:
# #                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
# #   ENC =~                                                                
# #     v5                0.995       NA                      0.962    0.814
# #     v7                0.850       NA                      0.822    0.768
# #     v17               1.110       NA                      1.074    0.886
# #   BUI =~                                                                
# #     v42r              0.864       NA                      0.777    0.672
# #     v43               0.689       NA                      0.620    0.571
# #   PSY =~                                                                
# #     v44               1.053       NA                      1.118    0.901
# #     v45               0.623       NA                      0.662    0.551
# # 
# # Regressions:
# #                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
# #   PSY ~                                                                 
# #     ENC              -0.110       NA                     -0.100   -0.100
# #     BUI               1.003       NA                      0.850    0.850
# # 
# # Covariances:
# #                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
# #   ENC ~~                                                                
# #     BUI              -0.264       NA                     -0.303   -0.303


```

  
#### Identification  
  
Before proceeding with the analysis, we need to check if our model is underidentified, just-identified or overidentified (Kline, 2016, pp.127-128). Model specifications can not require more parameters to be analysed than there are observations available. Thus, only *just-identified* (equal number of observations and parameters) and *overidentified* (less parameters than observations) models can be analysed (underidentification produces an error message).  
  
The number of observations (free parameters) to be estimated (non-reduntant elements) can be calculated with k(k+1)/2. "k" is the number of observed (manifest) variables in the model.      
  
Figure 4 above shows that we have seven observed manifest variables in our model (v5, v7, v17, v42r, v43, v44, v45). If we place that number into "k" and calculate, we get 7(7+1)/2 = 56/2 = 28. So, our number of observations is 28.  Next, we need to figure out how many parameters are to be estimated in the model (Figure 4). We will use *lavaan::lavInspect()* to analyse this:  
  
  
```{r number of parameters in the model, warning = F, message = F}

# how many parameters are to be estimated in the model?
# lavaan::inspect(sem_model_fit) # n of parameters: 24 # older version
lavaan::lavInspect(sem_model_fit) # n of parameters: 24

```
  
Result shows that there are 24 parameters to be estimated:  
  
4 manifest variable loadings on three factors (lambda matrix)  
7 error variances of manifest variables (theta matrix)  
3 latent variances and 1 covariance (psi matrix)  
2 regressions (beta matrix)  
7 intercepts(nu matrix)  
  
The model is overidentified (good!), as the number of non-redundant elements [7(7+1)/2 = 56/2 = 28] is greater than the number of parameters to be estimated: 28 - 24 = 4.  
  
  
#### Sample size  
  
Another issue is that do we have enough data for the SEM model? According to Jackson (2003), we can get a rough estimate of the sample size needed by multiplying the number of participants (N) by the number of parameters in the model (q). Sample-size-to-parameters ratio according to this ***N:q*** rule would ideally be 20:1 (Kline, 2016, p.16).  
  
We have 24 parameters to be estimated, so our sample size should be 20*24 = 480. Our data frame has `r dec0(nrow(dfttu))` rows (=participants), so we do fine here. However, Kline (2016, p.16) mentions that "Less ideal would be an N:q ratio of 10:1 ...", which in our case would mean radically smaller sample size expectation (of 240).  
  
A more clever way of investigating sample size needed is to use simulation alongside with desired alpha and power values. Usually alpha is set to 0.05 (five percent chance of Type I error) and power to 0.80 (80 percent chance to avoid rejecting H1 when it is true).  
  
```{r simple sample size detection}

# install.packages("semPower")
library(semPower)

aPrioriSampleSize <- semPower::semPower.aPriori(effect = .05, effect.measure = 'RMSEA', alpha = .05, power = .80, df = 24)
summary(aPrioriSampleSize)
 # F0                        0.060000
 # RMSEA                     0.050000
 # Mc                        0.970446
 #                                   
 # df                        24      
 # Required Num Observations 375     
 #                                   
 # Critical Chi-Square       36.41502
 # NCP                       22.44000
 # Alpha                     0.050000
 # Beta                      0.201086
 # Power (1-beta)            0.798914
 # Implied Alpha/Beta Ratio  0.248650

```
  
We learn that if we estimate 24 parameters (alpha = .05, power = .80), required sample size would be 375.      
  
Let's take another opinion from *simsem* ...  
  
```{r simsem}

# install.packages("simsem")
library(simsem)
# https://github.com/simsem/simsem/wiki/Example-18:-Simulation-with-Varying-Sample-Size

# define population model
sem_population_model <- '

# Factors
ENC =~ 0.81*v5   + 0.77*v7 + 0.89*v17   # standardized loadings from "Std.all" output above
BUI =~ 0.67*v42r + 0.57*v43
PSY =~ 0.90*v44  + 0.55*v45

# Regressions
PSY ~ -0.10*ENC + 0.85*BUI              # standardized path coefficients from "Std.all" output above

# Covariances
ENC ~~ -.30*BUI                         # covariance from "Std.all" output above
'

# simulation takes some time, so we'll save the results into a file ...
if (file.exists("SEM_sampleSize.rds")) {
SEM_sampleSize <- readRDS("SEM_sampleSize.rds")
} else {
SEM_sampleSize <- simsem::sim(NULL, n=50:500, model = sem_model, generate = sem_population_model, std.lv = TRUE, lavaanfun = "cfa")
saveRDS(SEM_sampleSize, "SEM_sampleSize.rds")
}

summary(SEM_sampleSize)

# look at alpha value at 0.05 and sample sizes
simsem::plotCutoff(SEM_sampleSize, 0.05)

simsem::getCutoff(SEM_sampleSize, 0.05, nVal = 240)	# 10:1 rule
#      chisq      aic      bic      rmsea       cfi       tli       srmr
# 1 19.40819 5487.851 5546.762 0.06378119 0.9449914 0.8949836 0.05273879

simsem::getCutoff(SEM_sampleSize, 0.05, nVal = 480)	# 20:1 rule
#      chisq      aic      bic      rmsea       cfi       tli       srmr
# 1 18.35418 10912.44 10996.96 0.03367612 0.9878986 0.9768973 0.02682591

```
  
```{r sample size needed for 0.80 power}

# find the sample size that provides the power just over 0.80
Cpow <- simsem::getPower(SEM_sampleSize)
simsem::findPower(Cpow, "N", 0.80)
   # ENC=~v5    ENC=~v7   ENC=~v17  BUI=~v42r   BUI=~v43   PSY=~v44   PSY=~v45    PSY~ENC    PSY~BUI   ENC~~BUI     v5~~v5     v7~~v7   v17~~v17 v42r~~v42r   v43~~v43   v44~~v44   v45~~v45 
   #     Inf        Inf        Inf         68         68        142        101         NA        246        266         50         54         65         92         50        306         58

```

  
```{r plot power against sample size}

simsem::plotPower(SEM_sampleSize, powerParam=c("ENC=~v5", "BUI=~v43", "PSY=~v44", "PSY~ENC", "PSY~BUI", "ENC~~BUI"))

```

  
Results of the simulation show that no sample size is enough to make "PSY~ENC" regression plausible if the power level is set to .80. However, investigation of "PSY~BUI" regression would require 246 participants (result in accordance with earlier "rule-of-thumb" 10:1 estimation). The greatest suggested sample size is 306, quite close to the earlier estimation of 375.    
  
We conclude that our model is overidentified and that we have enough participants.  
  
  
  
### 7.3 Extract the results  

```{r produce sem model summary, message = FALSE, warning = FALSE}

summary(sem_model_fit, standardized=TRUE, rsq=TRUE) # display R Squared

```
  
  
```{r optional results output, message = FALSE, warning = FALSE}

# https://rpubs.com/mkearney/103040

library(stargazer)

# unstandardised
pars.factors <- parameterEstimates(sem_model_fit)[ parameterEstimates(sem_model_fit)[,'op']=='=~', c(1:5)]
stargazer(pars.factors, summary=FALSE, type='text', # use 'html' if you wish to save the table and edit it further in Word etc.
          rownames=FALSE, initial.zero=FALSE, digits=3, title='Factor loadings')

pars.regressions <- parameterEstimates(sem_model_fit)[ parameterEstimates(sem_model_fit)[,'op']=='~', c(3,4,5,7)]
stargazer(pars.regressions, summary=FALSE, type='text', rownames=FALSE, initial.zero=FALSE, digits=3, title='Predicting PSY')

# standardised
std.pars.factors <- standardizedSolution(sem_model_fit)[ standardizedSolution(sem_model_fit)[,'op']=='=~', c(1:5)]
stargazer(std.pars.factors, summary=FALSE, type='text', rownames=FALSE, initial.zero=FALSE, digits=3, title='Factor loadings (std)')

std.pars.regressions <- standardizedSolution(sem_model_fit)[ standardizedSolution(sem_model_fit)[,'op']=='~', c(3,4,5,7)]
stargazer(std.pars.regressions, summary=FALSE, type='text', rownames=FALSE, initial.zero=FALSE, digits=3, title='Predicting PSY (std)')

# produce joint tables
# display both un/standardised parameter estimates in a same table
joint.pars.factors <- cbind(pars.factors, std.pars.factors[4:5]) # pick only last two columns from std parameter estimates
names(joint.pars.factors)[7] <- "se.std" # rename the second "se" in column 7 as "se.std"
stargazer(joint.pars.factors, summary=FALSE, type='text', rownames=FALSE, initial.zero=FALSE, digits=3, title='Factor loadings')
# display both un/standardised regressions in a same table
joint.pars.regressions <- cbind(pars.regressions, std.pars.regressions[2:4]) # pick three columns from std regressions
names(joint.pars.regressions)[6] <- "se.std" # rename the second "se" in column 6 as "se.std"
names(joint.pars.regressions)[7] <- "pvalue.std" # rename the second "pvalue" in column 7 as "pvalue.std"
stargazer(joint.pars.regressions, summary=FALSE, type='text', rownames=FALSE, initial.zero=FALSE, digits=3, title='Predicting PSY')

```

  
  
  
### 7.4 Show fit measures  
  
*Lavaan* produces a whole bunch of fit measures:  
  
```{r all fit measures, message = FALSE, warning = FALSE}

fitMeasures(sem_model_fit)

```
  
```{r fit measures for inline text, message = FALSE, warning = FALSE}

# produce fit indices for inline text
sem.fits <- as.numeric(fitMeasures(sem_model_fit, c("chisq",            # sem.fits[1]
                                                    "df",               # sem.fits[2]
                                                    "pvalue",           # sem.fits[3] 
                                                    "cfi",              # sem.fits[4] 
                                                    "tli",              # sem.fits[5] 
                                                    "aic",              # sem.fits[6]
                                                    "bic",              # sem.fits[7]
                                                    "rmsea",            # sem.fits[8]
                                                    "rmsea.ci.lower",   # sem.fits[9]
                                                    "rmsea.ci.upper",   # sem.fits[10]
                                                    "rmsea.pvalue",     # sem.fits[11]
                                                    "srmr",             # sem.fits[12]
                                                    "baseline.chisq",   # sem.fits[13]
                                                    "baseline.df"       # sem.fits[14]
                                                    )))

```

```{r display fit measures, message = FALSE, warning = FALSE}

# create vector of desired measure names
fit.names <- c("chisq", "df", "pvalue", "cfi", "tli", "aic", "bic", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper", "rmsea.pvalue", "srmr", "baseline.chisq", "baseline.df")  
# merge names and values (from the previous "sem.fits") into df
df.fit.measures <- data.frame(fit.names,sem.fits)
names(df.fit.measures)[1] <- "Measure"
names(df.fit.measures)[2] <- "Value"
# use Stargazer to produce a table
stargazer(df.fit.measures, summary=FALSE, type='text', # use 'html' if you wish to save the table and edit it further in Word etc.
          rownames=FALSE, initial.zero=FALSE, digits=3, title='Fit measures')

```
  
  
Here are short descriptions of selected fit measures:  
  
**chisq** (Chi square).  
  
**df** (degrees of freedom) indicate free (not fixed) parameters of the model.  
  
**pvalue** (probability value) is the probability of getting as large a discrepancy as occurred with the present sample (under appropriate distributional assumptions and assuming a correctly specified model) (Arbuckle & Wothke, 1999). It tests the hypothesis that the model fits perfectly to the population. The value should be greater than pre-set probability level (usually .05). If the value is .05 or less, the departure of the data from the model is significant at the .05 level (model does not fit into data). However, the appropriateness of hypothesis testing in model fitting, even when the necessary distributional assumptions are met, is routinely questioned.  
  
**NFI** (Normed Fit Index) tells how big discrepancy there is between the model being evaluated (default model) and the baseline model (terribly fitting ‘independence model’). **TLI, CFI, RNI** (Tucker-Lewis Index, Comparative fit index, Relative Noncentrality Index) similarly compare the proposed model to a baseline model that all other models should be expected to exceed (Hair, Anderson, Tatham & Black, 1995, p. 685). Values close to one indicate a very good fit. According to Bentler & Bonett (1980, p. 600), referring to both the NFI and the TLI, “... models with overall fit indices of less than .90 can usually be improved substantially.”    
  
**RMSEA** (Root Mean Square Error of Approximation). The RMSEA is designed to evaluate the approximate fit of the model in the population (Kaplan, 2000, p.112). This indice is getting smaller as the df increases. In practice this means that models with large RMSEA values (e.g., 0.12) simplify the reality. This error could be estimated as follows (Browne & Cudeck, 1993; Kaplan, 2000, p. 113): < 0.05 ‘close fit’, 0.05 – 0.08 ‘fair fit’, 0.081 – 0.10 ‘mediocre fit’, > 0.10 ‘poor fit’.  
  
**SRMR** (Standardized Root Mean square Residual). The SRMR is a measure of the mean absolute correlation residual, the overall difference between the observed and predicted correlations (Kline, 2016, p. 277). This "badness-of-fit" statistic (Kline, 2016) indicates if the square root of the discrepancy between the sample and model covariance matrixes is too big.  [Cangur and Ercan (2015)](https://digitalcommons.wayne.edu/cgi/viewcontent.cgi?article=1810&context=jmasm) have collected synthesis of several fit measures: An acceptable SRMR value should be below .10, while a good value is less than .05.   
  
  
Next, we list the minimum set of fit statistics (Kline, 2016, p. 269) that should be reported alongside with the model results.  
  
### 7.4.1 Measures of absolute fit

*Absolute fit measures* tell how well the model explains the data (Kline, 2016, p. 266). The do that by determining the degree to which the model predicts the observed correlation matrix (Hair et al., 1995, p. 683).   
  
First, relative $\chi^2$ is calculated by dividing the $\chi^2$ with *df*, resulting in `r dec2(sem.fits[1]/sem.fits[2])`. Usually, values less than five are considered to be adequate (Marsh & Hocevar, 1985). However, some researchers argue that the value should be less than two (Byrne, 1989).  
  
Second, the RMSEA estimate of `r dec3(sem.fits[8])` is in the upper bound of the mediocre fit level (.08 - .10), indicating that the model over simplifies the reality. The upper limit of the 90 per cent confidence interval (`r dec3(sem.fits[10])`) supports the assumption that we have quite poor fitting model here.  
  
Third, SRMR estimate of `r dec3(sem.fits[12])` indicates that the discrepancy between the sample and model covariance matrix is not too big (as the value is below .10).  
  
All three absolute fit measures are singing the same song: The model is not poor, but it is (perhaps due to being so simple) not doing best justice to the sample.  
  
  
### 7.4.2 Measures of incremental (or relative) fit  
  
*Incremental fit measures* compare the proposed model (default model) to a baseline model (independence model) that all other models should be expected to exceed. Here we discuss only one measure, Comparative fit index (CFI).   
  
The CFI value of `r dec3(sem.fits[4])` is calculated by dividing the relative $\chi^2$ of the proposed model by the relative $\chi^2$ of the baseline model and subtracting the output from 1:  
  
CFI = 1 - ($\chi^2$<sub>proposed</sub> - *df*<sub>proposed</sub>) / ($\chi^2$<sub>baseline</sub> - *df*<sub>baseline</sub>)  
    = 1 - (chisq-df)/(baseline.chisq-baseline.df)  
    = `r dec3(1-(sem.fits[1]-sem.fits[2])/(sem.fits[13]-sem.fits[14]))`  
  
CFI value shows that the proposed model has a discrepancy that is `r dec2(sem.fits[4]*100)` per cent of the way between the (terribly fitting) independence model (CFI = 0) and the (perfectly fitting) saturated model (CFI = 1).   
  
Quite high incremental fit measure (CFI) value of `r dec2(sem.fits[4])` indicates that this model has clearly better fit than the poor fitting baseline (independence) model.  
  
  
  
### 7.5 Report the results  
  
```{r extract sem model estimates for inline text, message = FALSE, warning = FALSE}

# we report only stat sig results in the text (PSY ~ BUI)
sem_m5_bui_psy_est <-  parameterEstimates(sem_model_fit)$est[9]
sem_m5_bui_psy_se <-  parameterEstimates(sem_model_fit)$se[9]
sem_m5_bui_psy_p <-  parameterEstimates(sem_model_fit)$pvalue[9]

# R Square
sem_m5_psy_r2 <- lavaan::inspect(sem_model_fit,'r2')[8] # use inspect(sem_model_fit,'r2') to see how many columns the R^2 table has (we are interested in the last value for the PSY)

```

  
Results showed that Build-up of work requirements was statistically significantly positively related to the Psychic stress of the work, *B* = `r dec3(sem_m5_bui_psy_est)`, *SE* = `r dec3(sem_m5_bui_psy_se)`, *p* = `r dec3(sem_m5_bui_psy_p)`. Encouraging leadership had, as expected, a negative relation to the Psychic stress of the work, but it was not significant. This model (see Figure 5) explained (almost solely by the Build-up of work requirements) `r dec2(sem_m5_psy_r2*100)`% of the variance of the Psychic stress of the work.  
  
Fit measures are presented in Table 1. The model explained the data quite well (RMSEA = `r dec3(sem.fits[8])`, *C.I.* = `r dec3(sem.fits[9])` - `r dec3(sem.fits[10])`; SRMR = `r dec3(sem.fits[12])`). However, the *p*-value of $\chi^2$ test was below the desired .05 value (*p* = `r dec3(sem.fits[3])`) and the relative $\chi^2$ (`r dec2(sem.fits[1]/sem.fits[2])`) was above the desired value of 5. Incremental fit measure (CFI = `r dec3(sem.fits[4])`) showed that the model is superior compared to a baseline model.   
  
  
Figure 5. Relationship between Build-up of work requirements (BUI), Encouraging leadership (ENC) and Psychic stress of the work (PSY) (unstandardized estimates). 
  
```{r unstandardized values diagram sem model, message = FALSE, warning = FALSE}

semPaths(sem_model_fit, 
         layout="tree",
         "par", # unstandardized parameter estimates
         edge.label.cex=0.7,
         edge.width=0.2,
         edge.color="black",
         optimizeLatRes=TRUE,
         label.scale=TRUE,
         normalize=TRUE,
         fade=FALSE,
         as.expression=c("nodes"))


```
  
Table 1. Fit measures  
  
```{r table 1 display fit measures, message = F, echo = F}

sem.fits.table1 <- as.numeric(fitMeasures(sem_model_fit, c("chisq",            # sem.fits.table1[1]
                                                           "df",               # sem.fits.table1[2]
                                                           "pvalue",           # sem.fits.table1[3] 
                                                           "rmsea",            # sem.fits.table1[4]
                                                           "rmsea.ci.lower",   # sem.fits.table1[5]
                                                           "rmsea.ci.upper",   # sem.fits.table1[6]
                                                           "srmr",             # sem.fits.table1[7]
                                                           "cfi"               # sem.fits.table1[8] 
                                                           )))

fit.names.table1 <- c("chisq", "df", "pvalue", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper", "srmr", "cfi")  
df.fit.measures.table1 <- data.frame(fit.names.table1,sem.fits.table1)
names(df.fit.measures.table1)[1] <- "Measure"
names(df.fit.measures.table1)[2] <- "Value"
stargazer(df.fit.measures.table1, summary=FALSE, type='text', # use 'html' if you wish to save the table and edit it further in Word etc.
          rownames=FALSE, initial.zero=FALSE, digits=3)

```
  
  
### 7.6 Can we make the model better?  
  
Modification indices (MI) demonstrate how re-specification of the model may improve the model fit [[read more](https://benwhalley.github.io/just-enough-r/modification-indices.html)]. However, any modification of the model (according to modification indices) should be in accordance with the theoretical assumptions behind the model (see critical discussion about MI's in Hermida, 2015). 
  
Table below shows the first 15 MI values in descending order.   
  
  
```{r modification indices, message = FALSE, warning = FALSE}

modindices(sem_model_fit, # or modificationindices(sem_model_fit) 
           sort = TRUE, 
           maximum.number = 15) # show the 15 highest MI values

```
  
The first row represents a modification (MI = `r dec3(modindices(sem_model_fit, sort = TRUE)$mi[1])`) that is expected to have the strongest impact on the $\chi^2$ value of the model. In this case, it is suggested that we add a covariance (~~) between v43 ("My working pace has increased in recent years") and v44 ("I feel that I am beginning to dislike my work"). As we recall from the earlier investigation of the correlation matrix, the correlation between these two variables was quite high (*r* = `r dec3(lavCor(dfttu4cor)[1,7])`).  
  
Although these two observed variables load on two different latent variables (PSY v44; BUI v43), they are quite close to each other in theoretical sense (as both PSY and BUI measure negative experiences towards work). If this covariance is added to the model specification, previous $\chi^2$ value of `r dec3(sem.fits.table1[1])` would improve approximately by `r dec3(modindices(sem_model_fit, sort = TRUE)$mi[1])`). On the other hand, as Hermida (2015) argues, one should be critical towards the initial model (and its psychometric properties) and ask "Why the association between these two observed variables was not made visible when specifying the original model?"  
  
For the demonstrational purposes we re-specify the original model according to the highest MI value.  
  
```{r modify the sem model, message = FALSE, warning = FALSE}

# specify the modified model
sem_model_mod <- '
# Factors
ENC =~ NA*v5 + v7 + 1*v17   # Encouraging leadership
BUI =~ 1*v42r + v43         # Build-up of work requirements
PSY =~ 1*v44 + v45          # Psychic stress of the work
# Regressions
PSY ~ ENC + BUI             # Predict PSY with ENC and BUI
# Covariances
ENC ~~ BUI                  # Let ENC and BUI correlate
v43 ~~ v44                  # Covariance added by modification indice recommendation
'

# estimate the modified model
sem_model_mod_fit <- sem(sem_model_mod, 
                    data=dfttu4,
                    group = NULL, # no grouping variable (e.g., gender, control/intervention,...)
                    cluster = NULL, # no clustering variable (e.g., multilevel data may contain several values for each participant)
                    missing = "fiml.x" # full information maximum likelihood approach
                    )

# produce diagram of the model
semPaths(sem_model_mod_fit,
         layout="tree",
         "std", # standardized parameter estimates
         edge.label.cex=0.7,
         edge.width=0.2,
         edge.color="black",
         optimizeLatRes=TRUE,
         label.scale=TRUE,
         normalize=TRUE,
         fade=FALSE,
         as.expression=c("nodes"))

```
  
Figure above shows that now there is a connecting arc between v43 and v44.  
  
Here are the estimates of the modified model:  
  
```{r modified sem model estimates, message = FALSE, warning = FALSE}

summary(sem_model_mod_fit, standardized=TRUE, rsq=TRUE) # display R Squared

```
  
We can see that the $\chi^2$ ("Test statistic") value has become smaller and the *df* has also dropped by one (as we asked to estimate one additional covariance). Interestingly, the overall explanation of the PSY variance (R Square) has after the modification dropped from the original model's (quite high) value of `r dec2(sem_m5_psy_r2)` down to `r dec2(lavaan::inspect(sem_model_mod_fit,'r2')[8])`.   
  
Here is the comparison of the regressions of the original and modified models:  
  
```{r, message = FALSE, warning = FALSE}

std.pars.regressions <- standardizedSolution(sem_model_fit)[ standardizedSolution(sem_model_fit)[,'op']=='~', c(3,4,5,7)]
std.pars.regressions.mod <- standardizedSolution(sem_model_mod_fit)[ standardizedSolution(sem_model_mod_fit)[,'op']=='~', c(3,4,5,7)]

joint.pars.regressions.comparison <- cbind(std.pars.regressions, std.pars.regressions.mod[2:4]) # pick three columns from mod regressions
names(joint.pars.regressions.comparison)[2] <- "est.orig"
names(joint.pars.regressions.comparison)[3] <- "se.orig"
names(joint.pars.regressions.comparison)[4] <- "p.orig"
names(joint.pars.regressions.comparison)[5] <- "est.mod"
names(joint.pars.regressions.comparison)[6] <- "se.mod"
names(joint.pars.regressions.comparison)[7] <- "p.mod"
stargazer(joint.pars.regressions.comparison, summary=FALSE, type='text', rownames=FALSE, initial.zero=FALSE, digits=3, title='Predicting PSY (std regressions)')

```
  
We see that the regression on ENC has become statistically significant (p.mod column) as we have included a covariance between variables v43 and v44.  
  
Here is the comparison of the fit indices of the original ("Value") and modified ("Value.mod") models:  
  
```{r compare fit indices of original and modified models, message = FALSE, warning = FALSE}

sem.fits.mod <- as.numeric(fitMeasures(sem_model_mod_fit, c("chisq",           # sem.fits.table1[1]
                                                           "df",               # sem.fits.table1[2]
                                                           "pvalue",           # sem.fits.table1[3] 
                                                           "rmsea",            # sem.fits.table1[4]
                                                           "rmsea.ci.lower",   # sem.fits.table1[5]
                                                           "rmsea.ci.upper",   # sem.fits.table1[6]
                                                           "srmr",             # sem.fits.table1[7]
                                                           "cfi"               # sem.fits.table1[8] 
                                                           )))

df.fit.measures.mod <- data.frame(fit.names.table1,sem.fits.mod)
names(df.fit.measures.mod)[2] <- "Value.mod"

joint.fit.comparison <- cbind(df.fit.measures.table1, df.fit.measures.mod[2]) 
stargazer(joint.fit.comparison, summary=FALSE, type='text', rownames=FALSE, initial.zero=FALSE, digits=3, title='Comparison of fit indices')

```
  
As the table above shows, the fit indices of the modified model are considerably better:  
  
* chisq/df ratio is now below the suggested value of 5 (`r dec2(sem.fits.mod[1])` / `r dec0(sem.fits.mod[2])` = `r dec1(sem.fits.mod[1]/sem.fits.mod[2])`)    
* RMSEA and SRMR are now smaller (as desired)  
* CFI is closer to desired value of 1   
  
## Exercise 4  
  
The data frame contains also the fourth factor (latent variable), "Know-how rewarding" (REW). It is measured with two variables:  
v13 The organization rewards its employees’ professional knowledge and skills.  
v14 Employees with increased knowledge are given extra responsibility.  
  
Add "Know-how rewarding" (REW) to the previous model as an exogenous variable (alongside with ENC and BUI). Produce standardized results in a similar manner than in section "7.5 Report the results". That is, write text describing findings and produce figure and table similar to Figure 5 and Table 1.   
  
  
### References  
  
  
Beaujean, A. A. (2014). *Latent Variable Modeling Using R. A Step-by-Step Guide.* Routledge.  
  
Bentler, P. M., & Bonett, D. G. (1980). Significance tests and goodness of fit in the analysis of covariance
structures. *Psychological Bulletin*, *88*, 588–606.  
  
Bollen, K. (1989). *Structural Equations with Latent Variables*. New York: John Wiley & Sons.  
  
Cangur, S., & Ercan, I. (2015). Comparison of Model Fit Indices Used in Structural Equation Modeling Under Multivariate
Normality. *Journal of Modern Applied Statistical Methods*, *14*(1), Article 14. https://digitalcommons.wayne.edu/jmasm/vol14/iss1/14/
  
Hair, J. F., Anderson, R. E., Tatham R. L., & Black, W. C. (1995). *Multivariate Data Analysis*. Fourth edition. Prentice Hall.  
  
Hermida, R. (2015). The problem of allowing correlated errors in structural equation modeling: concerns and considerations. *Computational Methods in Social Sciences*, *3*(1), 5-17.
  
Jöreskog, K. G. (1969). A general approach to confirmatory maximum likelihood factor analysis. *Psychometrika*, *34*, 183–202.  
  
Kaplan, D. (2000). *Structural Equation Modeling*. Thousand Oaks: Sage.  
  
Kline, R. B. (2016). *Principles and practice of structural equation modeling*. Fourth edition. Guilford Press. 
  
Marsh, H. W., & Hocevar, D. (1985). Application of confirmatory factor analysis to the study of self-concept: First- and higher-order factor models and their invariance across groups. *Psychological Bulletin*, *97*, 562–582.  
  
Nokelainen, P., & Ruohotie, P. (2009). Non-linear Modeling of Growth Prerequisites in a Finnish Polytechnic Institution of Higher Education. *Journal of Workplace Learning*, *21*(1), 36-57.  
  
Raykov, T., & Marcoulides, G. A. (2000). *A First Course in Structural Equation Modeling*. Lawrence Erlbaum Associates.  
  
Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling. *Journal of Statistical Software*, *48*(2), 1-36.  [http://www.jstatsoft.org/v48/i02/](http://www.jstatsoft.org/v48/i02/)  
  
Schumacker, R. E., & Lomax, R. G. (2004). *A Beginner’s Guide to Structural Equation Modeling*. Second edition. Lawrence Erlbaum Associates.  
  
Suhr, D. (2006). [The Basics of Structural Equation Modeling]((https://www.lexjansen.com/wuss/2006/tutorials/TUT-Suhr.pdf).  
  