---
title: "Controlling for Variables"
author: "Donny Williams"
date: "5/25/2020"
bibliography: ../inst/REFERENCES.bib
output:
rmarkdown::html_vignette:
toc: yes
vignette: >
%\VignetteIndexEntry{Controlling for Variables}
%\VignetteEngine{knitr::rmarkdown}
\usepackage[utf8]{inputenc}
---
# Introduction
This vignette describes how to control for variables. This is a new feature to **BGGM** (version `2.0.0`).
# Example 1: Multivariate Regression
When controlling for variables, a multivariate regression is fitted in **BGGM**. In fact, a GGM can be understood as a multivariate regression with intercepts only models for the predictors.
## Notes about Implementation
**BGGM** does not use the typical approach for multivariate regression in `R`. This avoids having to
write out each outcome variable, of which there are typically many in a GGM. In **BGGM**, it is assumed
that the data matrix includes only the variables to be included in the GGM and the control variables.
### Correct
Suppose that we want to control for education level, with five variables included in the
GGM.
```r
# data
Y <- bfi[,c(1:5, 27)]
# head
head(Y)
#> A1 A2 A3 A4 A5 education
#> 61617 2 4 3 4 4 NA
#> 61618 2 4 5 2 5 NA
#> 61620 5 4 5 4 4 NA
#> 61621 4 4 6 5 5 NA
#> 61622 2 3 3 4 5 NA
#> 61623 6 6 5 6 5 3
```
Notice that `Y` includes **only** the five variables and `education`.
## Fit Model
This model can then be fitted with
```
fit <- explore(Y, formula = ~ as.factor(education))
```
To show this is indeed a multivariate regression, here are the summarized regression coefficients for the first
outcome.
```
summ_coef <- regression_summary(fit)
# outcome one
summ_coef$reg_summary[[1]]
#> Post.mean Post.sd Cred.lb Cred.ub
#> (Intercept) 0.256 0.095 0.072 0.442
#> as.factor(education)2 0.073 0.128 -0.177 0.323
#> as.factor(education)3 -0.202 0.104 -0.405 -0.001
#> as.factor(education)4 -0.462 0.119 -0.691 -0.233
#> as.factor(education)5 -0.578 0.117 -0.815 -0.346
```
And here are the coefficients from `lm` (a univariate regression for `A1`)
```
round(
cbind(
# summary: coef and se
summary( lm(scale(A1, scale = F) ~ as.factor(education), data = Y))$coefficients[,1:2],
# confidence interval
confint( lm(scale(A1, scale = F) ~ as.factor(education), data = Y))
), 3)
#> Estimate Std. Error 2.5 % 97.5 %
#> (Intercept) 0.256 0.093 0.073 0.438
#> as.factor(education)2 0.072 0.125 -0.172 0.316
#> as.factor(education)3 -0.203 0.101 -0.401 -0.004
#> as.factor(education)4 -0.461 0.116 -0.690 -0.233
#> as.factor(education)5 -0.578 0.115 -0.804 -0.351
```
The estimate are very (very) similar.
## Summary
Note that all the other functions work just the same. For example, the relations controlling for education
are summarized with
```
summary(fit)
#> BGGM: Bayesian Gaussian Graphical Models
#> ---
#> Type: continuous
#> Analytic: FALSE
#> Formula: ~ as.factor(education)
#> Posterior Samples: 5000
#> Observations (n):
#> Nodes (p): 5
#> Relations: 10
#> ---
#> Call:
#> estimate(Y = Y, formula = ~as.factor(education))
#> ---
#> Estimates:
#> Relation Post.mean Post.sd Cred.lb Cred.ub
#> A1--A2 -0.239 0.020 -0.278 -0.200
#> A1--A3 -0.109 0.020 -0.150 -0.070
#> A2--A3 0.276 0.019 0.239 0.312
#> A1--A4 -0.013 0.021 -0.055 0.026
#> A2--A4 0.156 0.020 0.117 0.196
#> A3--A4 0.173 0.020 0.134 0.214
#> A1--A5 -0.010 0.020 -0.050 0.029
#> A2--A5 0.150 0.020 0.111 0.189
#> A3--A5 0.358 0.018 0.322 0.392
#> A4--A5 0.121 0.020 0.082 0.159
#> ---
```
### Incorrect
Now if we wanted to control for education, but also had gender in `Y`, this would be incorrect
```
Y <- bfi[,c(1:5, 26:27)]
head(Y)
#> A1 A2 A3 A4 A5 gender education
#> 61617 2 4 3 4 4 1 NA
#> 61618 2 4 5 2 5 2 NA
#> 61620 5 4 5 4 4 2 NA
#> 61621 4 4 6 5 5 2 NA
#> 61622 2 3 3 4 5 1 NA
#> 61623 6 6 5 6 5 2 3
```
In this case, with `estimate(Y, formula = as.factor(education))`, the GGM would also include `gender`
(six variables instead of the desired 5). This is because all variables not included in `formula` are included in the GGM. This was adopted in **BGGM** to save the user from having to write out each outcome.
This differs from `lm`, where each outcome needs to be written out, for example `cbind(A1, A2, A3, A4, A4) ~ as.factor(education)`. This is quite cumbersome for a model that includes many nodes.
# Example 2: Multivariate Probit
The above data is ordinal. In this case, it is possible to fit a multivariate probit model. This is also the approach for binary data in **BGGM**. This is implemented with
```
fit <- estimate(Y, formula = ~ as.factor(education),
type = "ordinal", iter = 1000)
```
Note that the multivariate probit models can also be summarized with `regression_summary`.
# Example 3: Gaussian Copula Graphical Model
This final example fits a Gaussian copula graphical model that can be used for mixed data. In this case,
`formula` is not used and instead all of the variables are included in the GGM.
## Fit Model
This model is estimated with
```
# data
Y <- na.omit(bfi[,c(1:5, 27)])
# fit type = "mixed"
fit <- estimate(Y, type = "mixed", iter = 1000)
# summary
summary(fit)
#> BGGM: Bayesian Gaussian Graphical Models
#> ---
#> Type: mixed
#> Analytic: FALSE
#> Formula:
#> Posterior Samples: 1000
#> Observations (n):
#> Nodes (p): 6
#> Relations: 15
#> ---
#> Call:
#> estimate(Y = Y, type = "mixed", iter = 1000)
#> ---
#> Estimates:
#> Relation Post.mean Post.sd Cred.lb Cred.ub
#> A1--A2 -0.217 0.048 -0.294 -0.114
#> A1--A3 -0.063 0.027 -0.113 -0.011
#> A2--A3 0.364 0.023 0.317 0.410
#> A1--A4 0.116 0.038 0.048 0.192
#> A2--A4 0.241 0.031 0.182 0.303
#> A3--A4 0.228 0.026 0.174 0.275
#> A1--A5 0.057 0.031 0.003 0.120
#> A2--A5 0.186 0.027 0.135 0.241
#> A3--A5 0.438 0.019 0.399 0.474
#> A4--A5 0.151 0.025 0.103 0.199
#> A1--education -0.016 0.069 -0.125 0.119
#> A2--education 0.063 0.049 -0.016 0.162
#> A3--education 0.049 0.025 0.002 0.099
#> A4--education 0.053 0.026 0.005 0.105
#> A5--education 0.072 0.024 0.024 0.120
#> ---
```
Here it is clear that education is included in the model, as the relations with the other nodes are included in the output.
## Select Graph
The graph is selected with
```
select(fit)
```
# Note
It is possible to control for variable with all methods in **BGGM**, including when comparing groups, Bayesian hypothesis testing, etc.