# Load Packages
library(lavaan)
library(gam)
library(Matrix)
library(expm)
# enter correlation matrix
rmatrix<-"
1
0.4 1
0.3 0.2 1
0.2 0.1 0.3 1
"
rmatrix<-getCov(rmatrix)
dimnames(rmatrix) <- list(c("IQ","Verbal","Maths","Space"), c("IQ","Verbal","Maths","Space"))
rmatrix
## IQ Verbal Maths Space
## IQ 1.0 0.4 0.3 0.2
## Verbal 0.4 1.0 0.2 0.1
## Maths 0.3 0.2 1.0 0.3
## Space 0.2 0.1 0.3 1.0
######### Extract the first principal component ###################
# Step 1: Sum the coefficients in each column
u1 <- matrix(1,1,4) %*% rmatrix
u1
## IQ Verbal Maths Space
## [1,] 1.9 1.7 1.8 1.6
# Step 2: Normalizing u1
t1 <- sqrt(u1 %*% t(u1))
t1
## [,1]
## [1,] 3.507136
# Computing the first trial vector
v1 <- u1/t1[1,1]
v1
## IQ Verbal Maths Space
## [1,] 0.5417527 0.4847261 0.5132394 0.4562128
# Step 3: Computing the second trial vector
u2 <- v1 %*% rmatrix
u2
## IQ Verbal Maths Space
## [1,] 0.9808574 0.8496963 0.9095742 0.7670077
# Step4: Normalizing u2
t2 <- sqrt(u2 %*% t(u2))
t2
## [,1]
## [1,] 1.760594
# Computing the second trial vector
v2 <- u2/t2[1,1]
v2
## IQ Verbal Maths Space
## [1,] 0.5571173 0.4826191 0.5166292 0.4356528
# Step 5: Comparing the first two trials
vectors1.2 <- rbind(v1,v2)
vectors1.2
## IQ Verbal Maths Space
## [1,] 0.5417527 0.4847261 0.5132394 0.4562128
## [2,] 0.5571173 0.4826191 0.5166292 0.4356528
# Step 6: Produce the third trial vector
u3 <- v2 %*% rmatrix
u3
## IQ Verbal Maths Space
## [1,] 0.9922843 0.8523572 0.910984 0.7503269
# Step 7: Normalizing u3
t3 <- sqrt(u3 %*% t(u3))
t3
## [,1]
## [1,] 1.761824
# Compute v3
v3 <- u3/t3[1,1]
v3
## IQ Verbal Maths Space
## [1,] 0.5632143 0.4837925 0.5170687 0.4258808
# Step 8: Compare the second and thirs trial vectors
vectors2.3 <- rbind(v2, v3)
vectors2.3
## IQ Verbal Maths Space
## [1,] 0.5571173 0.4826191 0.5166292 0.4356528
## [2,] 0.5632143 0.4837925 0.5170687 0.4258808
# Step 9: Compute the component loading for the first component
eigen1 <- t3[1,1] # First eigenvalue
sqrt.eigen1 <- sqrt(eigen1)
pc1 <- sqrt.eigen1*v3
pc1
## IQ Verbal Maths Space
## [1,] 0.7475752 0.6421557 0.6863245 0.5652874
# Percent of variance explained by pc1
eigen1/4
## [1] 0.440456
########### Extract the second principal component ##############
# Step 10: Compute the first residual correlation matrix
reproduced.r <- t(pc1) %*% pc1
reproduced.r
## IQ Verbal Maths Space
## IQ 0.5588687 0.4800597 0.5130792 0.4225948
## Verbal 0.4800597 0.4123640 0.4407272 0.3630025
## Maths 0.5130792 0.4407272 0.4710414 0.3879706
## Space 0.4225948 0.3630025 0.3879706 0.3195498
R1 <- rmatrix - reproduced.r
R1
## IQ Verbal Maths Space
## IQ 0.44113134 -0.08005969 -0.2130792 -0.2225948
## Verbal -0.08005969 0.58763601 -0.2407272 -0.2630025
## Maths -0.21307919 -0.24072724 0.5289586 -0.0879706
## Space -0.22259483 -0.26300254 -0.0879706 0.6804502
# Reflecting residual matrix
for(i in 3:4){
R1[i,] <- R1[i,] * (-1)
R1[,i] <- R1[,i] * (-1)
}
R1
## IQ Verbal Maths Space
## IQ 0.44113134 -0.08005969 0.2130792 0.2225948
## Verbal -0.08005969 0.58763601 0.2407272 0.2630025
## Maths 0.21307919 0.24072724 0.5289586 -0.0879706
## Space 0.22259483 0.26300254 -0.0879706 0.6804502
# Step 11: Sum the coefficients in each column of R1
ur1 <- matrix(1,1,4) %*% R1
ur1
## IQ Verbal Maths Space
## [1,] 0.7967457 1.011306 0.8947945 1.078077
# Step 12: Normalizing ur1
tr1 <- sqrt(ur1 %*% t(ur1))
tr1
## [,1]
## [1,] 1.902748
# Computing the first trial vector
vr1 <- ur1/tr1[1,1]
vr1
## IQ Verbal Maths Space
## [1,] 0.4187342 0.5314976 0.4702643 0.5665894
# Step 13: Computing the second trial vector
ur2 <- vr1 %*% R1
ur2
## IQ Verbal Maths Space
## [1,] 0.3684886 0.5410233 0.4160766 0.5771597
# Step 14: Normalizing u2
tr2 <- sqrt(ur2 %*% t(ur2))
tr2
## [,1]
## [1,] 0.9668108
# Computing the second trial vector
vr2 <- ur2/tr2[1,1]
vr2
## IQ Verbal Maths Space
## [1,] 0.3811383 0.5595958 0.4303599 0.5969727
# Step 15: Comparing the first two trials
vectors.r.1.2 <- rbind(vr1,vr2)
vectors.r.1.2
## IQ Verbal Maths Space
## [1,] 0.4187342 0.5314976 0.4702643 0.5665894
## [2,] 0.3811383 0.5595958 0.4303599 0.5969727
# Step 16: Produce the third trial vector
ur3 <- vr2 %*% R1
ur3
## IQ Verbal Maths Space
## [1,] 0.3479148 0.5589295 0.3910491 0.6003657
# Step 17: Normalizing u3
tr3 <- sqrt(ur3 %*% t(ur3))
tr3
## [,1]
## [1,] 0.9730392
# Compute vr3
vr3 <- ur3/tr3[1,1]
vr3
## IQ Verbal Maths Space
## [1,] 0.3575547 0.5744162 0.4018842 0.6170005
# Step 17: Produce the fourth trial vector
ur4 <- vr3 %*% R1
ur4
## IQ Verbal Maths Space
## [1,] 0.3347153 0.5679391 0.3727673 0.6151469
# Step 18: Normalizing ur4
tr4 <- sqrt(ur4 %*% t(ur4))
tr4
## [,1]
## [1,] 0.9756794
# Compute vr4
vr4 <- ur4/tr4[1,1]
vr4
## IQ Verbal Maths Space
## [1,] 0.3430587 0.582096 0.3820592 0.6304805
# Step 19: Compare the second and thirs trial vectors
vectors.r.3.4 <- rbind(vr3, vr4)
vectors.r.3.4
## IQ Verbal Maths Space
## [1,] 0.3575547 0.5744162 0.4018842 0.6170005
## [2,] 0.3430587 0.5820960 0.3820592 0.6304805
# Step 20: Produce the fifth trial vector
ur5 <- vr4 %*% R1
ur5
## IQ Verbal Maths Space
## [1,] 0.3264821 0.5723855 0.3598548 0.6248564
# Step 21: Normalizing ur4
tr5 <- sqrt(ur5 %*% t(ur5))
tr5
## [,1]
## [1,] 0.9768094
# Compute vr5
vr5 <- ur5/tr5[1,1]
vr5
## IQ Verbal Maths Space
## [1,] 0.3342331 0.5859745 0.3683982 0.6396912
# Step 22: Compare the fourth and fifth trial vectors
vectors.r.4.5 <- rbind(vr4, vr5)
vectors.r.4.5
## IQ Verbal Maths Space
## [1,] 0.3430587 0.5820960 0.3820592 0.6304805
## [2,] 0.3342331 0.5859745 0.3683982 0.6396912
# Step 23: Compute the component loading for the second component
eigen2 <- tr5[1,1] # First eigenvalue
sqrt.eigen2 <- sqrt(eigen2)
pc2 <- sqrt.eigen2*vr5
pc2
## IQ Verbal Maths Space
## [1,] 0.3303349 0.5791401 0.3641014 0.6322303
# Percent of variance explained by pc1
eigen2/4
## [1] 0.2442024
# Component loadings matrix
pc1 <- t(pc1)
pc2 <- t(pc2)
loading <- data.frame(pc1, pc2)
loading
## pc1 pc2
## IQ 0.7475752 0.3303349
## Verbal 0.6421557 0.5791401
## Maths 0.6863245 0.3641014
## Space 0.5652874 0.6322303