PS15 Solutions
Sarah
4/19/2022
STAT 360: Computational Statistics and Data Analysis
Load R Libraries, Import and Attach Relevant Data, and Specify Seed
library(rmarkdown); library(knitr); library(readxl)
set.seed(42069)EXERCISE 01
mentalCov <- matrix(c(.165, .031, .055, .038, .031, .205, 0.089, 0.040, .055, 0.089, .165, .053, .038, 0.040, .053, .172), nrow = 4, ncol = 4)
rownames(mentalCov)<- c("Comfortable", "Benefits", "Anonymity", "Discussed")
colnames(mentalCov)<- c("Comfortable", "Benefits", "Anonymity", "Discussed")
R <- solve(sqrt(diag(diag(mentalCov)))) %*% mentalCov %*%
t(solve(sqrt(diag(diag(mentalCov)))))
rownames(R)<- c("Comfortable", "Benefits", "Anonymity", "Discussed")
colnames(R)<- c("Comfortable", "Benefits", "Anonymity", "Discussed")
R## Comfortable Benefits Anonymity Discussed
## Comfortable 1.0000000 0.1685554 0.3333333 0.2255680
## Benefits 0.1685554 1.0000000 0.4839173 0.2130192
## Anonymity 0.3333333 0.4839173 1.0000000 0.3146079
## Discussed 0.2255680 0.2130192 0.3146079 1.0000000library(knitr)
include_graphics("C:/Users/Sarah Chock/OneDrive - University of St. Thomas/Senior Year/STAT 360 Comp Stat and Data Analysis/Structural Equation Models/PS14 Q1.png")
Part (a)
Independent variables: Discussed, Healthcare Plan, E1, E2, E3
Dependent Variables: Benefits, Anonymity, Comfort
PHI:phi <- diag(rep(.2,5))
rownames(phi)<- c("Discussed", "Healthcare Plan", "E1", "E2", "E3")
colnames(phi)<- c("Discussed", "Healthcare Plan", "E1", "E2", "E3")
phi## Discussed Healthcare Plan E1 E2 E3
## Discussed 0.2 0.0 0.0 0.0 0.0
## Healthcare Plan 0.0 0.2 0.0 0.0 0.0
## E1 0.0 0.0 0.2 0.0 0.0
## E2 0.0 0.0 0.0 0.2 0.0
## E3 0.0 0.0 0.0 0.0 0.2Part (b)
Gx <- matrix(0, nrow = 1, ncol = 5)
rownames(Gx)<- c("Discussed")
colnames(Gx)<- c("Discussed", "Healthcare Plan", "E1", "E2", "E3")
Gx[1,1] = 1
Gx## Discussed Healthcare Plan E1 E2 E3
## Discussed 1 0 0 0 0Part (c)
SigmaXX <- Gx %*% phi %*% t(Gx)
SigmaXX## Discussed
## Discussed 0.2Part (d)
So, we were hoping to get a value of .172, and we got .2. I think that isn't too shabby, but at this point we don't have a way to quantify how good or bad it is. I think it is fine. EXERCISE 02
Part (a)
beta <- matrix(0, nrow = 6, ncol = 6)
rownames(beta)<- c("Threats", "Assaults", "Fights", "Depression", "Anxiety", "Psychological Distress")
colnames(beta)<- c("Threats", "Assaults", "Fights", "Depression", "Anxiety", "Psychological Distress")
beta["Depression","Psychological Distress"] = 1
beta["Anxiety","Psychological Distress"] = .3
beta## Threats Assaults Fights Depression Anxiety
## Threats 0 0 0 0 0
## Assaults 0 0 0 0 0
## Fights 0 0 0 0 0
## Depression 0 0 0 0 0
## Anxiety 0 0 0 0 0
## Psychological Distress 0 0 0 0 0
## Psychological Distress
## Threats 0.0
## Assaults 0.0
## Fights 0.0
## Depression 1.0
## Anxiety 0.3
## Psychological Distress 0.0identity <- diag(c(1,1,1,1,1,1))
identity## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0 0 0 0
## [2,] 0 1 0 0 0 0
## [3,] 0 0 1 0 0 0
## [4,] 0 0 0 1 0 0
## [5,] 0 0 0 0 1 0
## [6,] 0 0 0 0 0 1gamma <- matrix(0, nrow = 6, ncol = 8)
rownames(gamma)<- c("Threats", "Assaults", "Fights", "Depression", "Anxiety", "Pyschological Distress")
colnames(gamma)<- c("Housing", "Violence", "E1", "E2", "E3", "E4", "E5", "D1")
gamma[6,1] = .7
gamma[6,2] = .7
gamma[6,8] = 1
gamma[5,7] = 1
gamma[4,6] = 1
gamma[3,2] = .7
gamma[3,5] = 1
gamma[2,2] = .7
gamma[2,4] = 1
gamma[1,2] = 1
gamma[1,3] = 1
gamma## Housing Violence E1 E2 E3 E4 E5 D1
## Threats 0.0 1.0 1 0 0 0 0 0
## Assaults 0.0 0.7 0 1 0 0 0 0
## Fights 0.0 0.7 0 0 1 0 0 0
## Depression 0.0 0.0 0 0 0 1 0 0
## Anxiety 0.0 0.0 0 0 0 0 1 0
## Pyschological Distress 0.7 0.7 0 0 0 0 0 1phi <- diag(rep(.5,8))
rownames(phi)<- c("Housing", "Violence", "E1", "E2", "E3", "E4", "E5", "D1")
colnames(phi)<- c("Housing", "Violence", "E1", "E2", "E3", "E4", "E5", "D1")
phi## Housing Violence E1 E2 E3 E4 E5 D1
## Housing 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## Violence 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0
## E1 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0
## E2 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0
## E3 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0
## E4 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0
## E5 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0
## D1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5Gy <- matrix(0, nrow = 5, ncol = 6)
rownames(Gy)<- c("Threats", "Assaults", "Fights", "Depression", "Anxiety")
colnames(Gy)<- c("Threats", "Assaults", "Fights", "Depression", "Anxiety", "Psychological Distress")
Gy[1,1] = 1
Gy[2,2] = 1
Gy[3,3] = 1
Gy[4,4] = 1
Gy[5,5] = 1
Gy## Threats Assaults Fights Depression Anxiety Psychological Distress
## Threats 1 0 0 0 0 0
## Assaults 0 1 0 0 0 0
## Fights 0 0 1 0 0 0
## Depression 0 0 0 1 0 0
## Anxiety 0 0 0 0 1 0Gx <- matrix(0, nrow = 1, ncol = 8)
rownames(Gx)<- c("Housing")
colnames(Gx)<- c("Housing", "Violence", "E1", "E2", "E3", "E4", "E5", "D1")
Gx[1,1] = 1
Gx## Housing Violence E1 E2 E3 E4 E5 D1
## Housing 1 0 0 0 0 0 0 0Part (b)
Sigmayy <- Gy %*% solve(identity - beta) %*% gamma %*% phi %*% t(gamma) %*% t(solve(identity - beta)) %*% t(Gy)
Sigmayy## Threats Assaults Fights Depression Anxiety
## Threats 1.000 0.3500 0.3500 0.350 0.1050
## Assaults 0.350 0.7450 0.2450 0.245 0.0735
## Fights 0.350 0.2450 0.7450 0.245 0.0735
## Depression 0.350 0.2450 0.2450 1.490 0.2970
## Anxiety 0.105 0.0735 0.0735 0.297 0.5891Part (c)
Sigmayx <- Gy %*% solve(identity - beta) %*% gamma %*% phi %*% t(Gx)
Sigmayx## Housing
## Threats 0.000
## Assaults 0.000
## Fights 0.000
## Depression 0.350
## Anxiety 0.105Part (d)
Sigmaxy <- t(Gy %*% solve(identity - beta) %*% gamma %*% phi %*% t(Gx))
Sigmaxy## Threats Assaults Fights Depression Anxiety
## Housing 0 0 0 0.35 0.105Part (e)
Sigmaxx <- Gx %*% phi %*% t(Gx)
Sigmaxx## Housing
## Housing 0.5Part (f)
yyandyx <- cbind(Sigmayy, Sigmayx)
xyandxx <- cbind(Sigmaxy, Sigmaxx)
sigma <- rbind(yyandyx, xyandxx)
sigma## Threats Assaults Fights Depression Anxiety Housing
## Threats 1.000 0.3500 0.3500 0.350 0.1050 0.000
## Assaults 0.350 0.7450 0.2450 0.245 0.0735 0.000
## Fights 0.350 0.2450 0.7450 0.245 0.0735 0.000
## Depression 0.350 0.2450 0.2450 1.490 0.2970 0.350
## Anxiety 0.105 0.0735 0.0735 0.297 0.5891 0.105
## Housing 0.000 0.0000 0.0000 0.350 0.1050 0.500Part (g)
library(cats)
set.seed(69)
here_kitty()
## meow