Chapter 8 Notes
8.1 APY
APY is application of blockwise matrix inversion, or wiki link
\[ \begin{bmatrix} P_{11} & P_{12}\\ P_{21} & P_{22} \end{bmatrix} ^{-1} = \begin{bmatrix} P_{11}^-1+P_{11}^{-1}P_{12}F^{-1}P_{21}P_{11}^{-1}&P_{11}^{-1}P_{12}F^{-1}\\ F^{-1}P_{21}P_{11}^{-1}&F^{-1} \end{bmatrix}\]
in which \(F=P_{22}-P_{21}P_{11}^{-1}P_{12}\) is assumed nonsingular.
A much original idea of APY resembles Quaas and Pollak (J Anim Sci, 1980, 51:1277-87), an example of which can be found in Lynch&Walsh page 759-61.
\(G^{-1}= \left[ \begin{matrix} I& -P_{cn}\\ 0&I \end{matrix} \right] \left[ \begin{matrix} G_{cc}^{-1}& 0\\ 0&M_{nn}^{-1} \end{matrix} \right] \left[ \begin{matrix} I& 0\\ -P_{nc}&I \end{matrix} \right] = \left[ \begin{matrix}G_{cc}^{-1}&-P_{cn}M_{nn}^{-1}\\ 0&M^{-1}_{nn} \end{matrix} \right] \left[ \begin{matrix}I&0\\ -P_{nc}&I \end{matrix} \right]\)
Citation: Misztal, I, 2016, Inexpensive computation of the inverse of the genomic relationship matrix in populations with small effective population size, \(Genetics\), 202:401-409.
in which \(G\) is the numerical relationship matrix, and \(P_{cn}=G_{cc}^{-1}G_{cn}\) and \(P_{nc}=G_{nc}G_{cc}^{-1}\). \(G_{cc}\), \(G_{cn}\), and \(G_{nc}\) (\(G_{cn}=G_{nc}\)) are, “core” to “core”, “core” to “non-core”, and “non-core” to “core” relationship matrix.
\(M_{nn}=diag(G_{nn})-diag(P_{cn}^TG_{cn})\)
8.1.1 A numerical example
This example is taken from appendix in Misztal’s paper.
\(G=\left[ \begin{matrix} 0.81 & 0 & 0 & 0.80 & -0.80\\ &0.81 & 0 & 0.80 &-0.80\\ & & 0.01& 0 & 0\\ & & & 1.61 & -1.60\\ symm. & & & & 1.61 \end{matrix} \right]\)
and \(G_{cc}^{-1}=\left[ \begin{matrix} 1.235 & 0\\ 0 & 1.235 \end{matrix} \right]\)
\(G_{cn}=G_{nc}^T=\left [ \begin{matrix} 0 & 0.80 & -0.80\\ 0 & 0.80 & -0.80 \end{matrix} \right]\)
and
\(P_{cn}=G_{cc}^{-1}G_{cn}=\left [ \begin{matrix} 0.00 & 0.988 & -0.988\\ 0.00 & 0.988 & -0.988 \end{matrix} \right ]\)
\(M_{nn}=diag(G_{nn})-diag(P_{cn}^TG_{cn}) = \left [ \begin{matrix}0.01 & &\\ & 1.61 & \\ & & 1.61 \end{matrix} \right] - \left [ \begin{matrix}0.00 & &\\ & 1.58 & \\ & & 1.58 \end{matrix} \right]=\left [ \begin{matrix}0.01 & &\\ & 0.03 & \\ & & 0.03 \end{matrix} \right]\)
and
\(G^{-1}=\left[ \begin{matrix} 66.8 & 65.5 & 0 & -33.1 & 33.1\\ &66.804 & 0 & -33.1 & 33.1\\ & & 100& 0 & 0\\ & & & 33.6 & 0\\ symm. & & & & 33.6 \end{matrix} \right]\)
This is APY inversed \(G\), and for a comparison, a regular inverse of \(G\) is quit different \(G^{-1}_{reg}=\left[ \begin{matrix} 40.6 & 39.4 & 0 & -19.9 & 19.9\\ &40.6 & 0 & -19.9 & 19.9\\ & & 100& 0 & 0\\ & & & 60.0 & 39.9\\ symm. & & & & 60.0 \end{matrix} \right]\)
GRM_C=matrix(c(0.81, 0, 0, 0.81), 2, 2, byrow = T)
GRM_N=matrix(0, 3, 3)
diag(GRM_N)=c(.01, 1.61, 1.61)
GRM_N[2,3]=GRM_N[3,2]=-1.6
GRM_CN=matrix(0, 2, 3)
GRM_CN[,2]=0.8
GRM_CN[,3]=-0.8
GRM=rbind(cbind(GRM_C, GRM_CN), cbind(t(GRM_CN), GRM_N))
print(GRM)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.81 0.00 0.00 0.80 -0.80
## [2,] 0.00 0.81 0.00 0.80 -0.80
## [3,] 0.00 0.00 0.01 0.00 0.00
## [4,] 0.80 0.80 0.00 1.61 -1.60
## [5,] -0.80 -0.80 0.00 -1.60 1.61## [,1] [,2] [,3] [,4] [,5]
## [1,] 40.64222 39.40765 0 -19.95012 19.95012
## [2,] 39.40765 40.64222 0 -19.95012 19.95012
## [3,] 0.00000 0.00000 100 0.00000 0.00000
## [4,] -19.95012 -19.95012 0 60.09975 39.90025
## [5,] 19.95012 19.95012 0 39.90025 60.09975IC=solve(GRM_C)
Pcn=IC %*% GRM_CN
Pnc=t(GRM_CN) %*% IC
Pcn_Gcn=t(Pcn) %*% GRM_CN
Mnn=diag(diag(GRM_N)-diag(Pcn_Gcn), nrow=nrow(GRM_N), ncol=nrow(GRM_N))
IMnn=solve(Mnn)
v1_1=rbind(IC, matrix(0, nrow=nrow(GRM_N), ncol=ncol(GRM_C)))
v1_2=rbind(-1*Pcn%*%IMnn, IMnn)
V1=cbind(v1_1, v1_2)
v2_1=rbind(diag(1, nrow=nrow(GRM_C), ncol=nrow(GRM_C)), -1*Pnc)
v2_2=rbind(matrix(0, nrow=nrow(GRM_C), ncol=nrow(GRM_N)), diag(1, nrow(GRM_N), ncol(GRM_N)))
V2=cbind(v2_1, v2_2)
IV=V1 %*% V2
print(IV) #apy## [,1] [,2] [,3] [,4] [,5]
## [1,] 66.80498 65.57041 0 -33.19502 33.19502
## [2,] 65.57041 66.80498 0 -33.19502 33.19502
## [3,] 0.00000 0.00000 100 0.00000 0.00000
## [4,] -33.19502 -33.19502 0 33.60996 0.00000
## [5,] 33.19502 33.19502 0 0.00000 33.60996## [,1] [,2] [,3] [,4] [,5]
## [1,] 8.100000e-01 9.572467e-16 0.00 0.800000 -0.800000
## [2,] 1.089201e-15 8.100000e-01 0.00 0.800000 -0.800000
## [3,] 0.000000e+00 0.000000e+00 0.01 0.000000 0.000000
## [4,] 8.000000e-01 8.000000e-01 0.00 1.610000 -1.580247
## [5,] -8.000000e-01 -8.000000e-01 0.00 -1.580247 1.6100008.2 LD score regression
source("~/R/MyLib/shotgun.R")
REP = 50
N = 500
M = 1000
h2 = 0.5
h2E = array(0, dim = c(2, REP))
for (rep in 1:REP) {
fq = runif(M, 0.05, 0.95)
Dl = array(0, dim = M)
Dl = runif(M, 0.8, 0.9)
Dl[seq(10, M, 10)] = 0
G = GenerateGenoDprime(fq, Dl[1:(length(Dl) - 1)], N)
b = rnorm(M, 0, sqrt(h2/M))
bv = G %*% b
Y = bv + rnorm(N, 0, 1)
chi = array(0, dim = M)
for (i in 1:M) {
sm = summary(lm(Y ~ G[, i]))
chi[i] = sm$coefficients[2, 3]^2
}
plot(chi, b^2 * fq * (1 - fq), pch = 16)
lds = array(0, dim = M)
for (i in 1:(M/10)) {
cg = cor(G[, ((i - 1) * 10 + 1):(i * 10)])
for (j in 1:10) {
lds[(i - 1) * 10 + j] = sum(cg[, j]^2)
}
}
ldSc = lm(chi ~ lds)
h2E[1, rep] = ldSc$coefficients[2] * M/N
h2E[2, rep] = var(bv)/var(Y)
}
plot(h2E[1, ], h2E[2, ])8.3 Matrix inversion wiki
8.3.1 2 X 2
The cofactor equation listed above yields the following result for 2 ?? 2 matrices. Inversion of these matrices can be done as follows
\[\text{A}^{-1}=\left[ \begin{matrix} a & b\\ c & d \end{matrix} \right] ^{-1} =\frac{1}{det (\text{A})}\left[ \begin{matrix} A & B\\ C & D \end{matrix} \right] =\frac{1}{ad-bc}\left[ \begin{matrix} d & -b\\ -c & a \end{matrix} \right] \] Where the scalar \(A\) is not to be confused with the matrix A.
8.3.2 3 X 3
\[\text{A}^{-1}=\left[ \begin{matrix} a & b & c\\ d & e & f\\ g & h & i \end{matrix} \right] ^{-1} =\frac{1}{det (\text{A})}\left[ \begin{matrix} A & B & C\\ D & E & F\\ G & H & I \end{matrix} \right] ^{T} =\frac{1}{det (\text{A})}\left[ \begin{matrix} A & D & G\\ B & E & H\\ C & F & I \end{matrix} \right] \]
in which \[\begin{matrix} A=(ei-fh), & D=-(bi-ch), & G = -(bf-ce), \\ B=-(di-fg), & E=(ai-cg), & H = -(af-cd), \\ C=(dh-eg), & F=-(ah-bg), & I = (ae-bd). \end{matrix} \]