Calculating partitioned covariance matrices using R

Let's say my data consists of 25 observations of four features, which are in groups of 2. So we'll call my variables $$x_1, x_2, y_1, y_2$$. We have a partitioned sample mean vector given by $$\begin{bmatrix}\bar{x}\\ \bar{y}\end{bmatrix}$$, where $$\bar{x} = \begin{bmatrix}\bar{x_1}\\ \bar{x_2}\end{bmatrix}$$, $$\bar{y} = \begin{bmatrix}\bar{y_1}\\ \bar{y_2}\end{bmatrix}$$. We can partition the covariance matrix $$S$$ by $$\begin{equation*} S = \begin{bmatrix}S_{xx}&S_{yx}\\S_{xy}&S_{yy}\end{bmatrix} \end{equation*}$$ where $$S_{xx}$$ is the covariance matrix of just the $$x$$ variables and $$S_{yy}$$ is the covariance matrix of just the $$y$$ variables. I'm trying to write an R script that will calculate $$S_{yx}$$ using just the input data. Clearly $$S_{yx}$$ consists of $$\begin{equation*} S_{yx} = \begin{bmatrix}Cov(x_1, y_1)&Cov(x_1, y_2)\\Cov(x_2, y_1)&Cov(x_2, y_2)\end{bmatrix} \end{equation*}$$ So if $$X$$ is the $$25 \times 2$$ matrix of $$x$$ data and $$Y$$ is the $$25 \times 2$$ matrix of $$Y$$ data, then my idea for calculating is is like so: $$\begin{equation*} S_{yx} = \frac{1}{n-1} \big(X^TY - n\bar{x}\bar{y}^T\big) \end{equation*}$$ However, when I type into R:

(1/(n-1)) * ((t(X) %*% Y) - (n * xbar %*% t(ybar)))


I get results which are wildly different from the corresponding results in the full covariance matrix. Is there a problem with my reasoning or my R code?

EDIT: ftp://ftp.wiley.com/public/sci_tech_med/multivariate_analysis_3e/ My input is the "T3_8_SONS.DAT" file in this directory. I don't know if there is a better way to include my input, sorry.

• It''s always helpful to include the inputs to your code. Feb 24, 2019 at 0:02
• @StatsStudent I updated my post with the input, thank you Feb 24, 2019 at 0:11
• I realized today too that I failed to actually address your question about the calculation you used to obtain $Syx.$ The formula you have is correct and indeed works, if $X$, $Y$, $\bar{X}$, and $\bar{Y}$ are properly sorted. I've provided an example in my answer. Feb 25, 2019 at 17:48

You can make your life a lot easier by using R's covariance function cov. So, all you need to do is read in the data to a dataframe, change the order of the columns, and then compute the covariance matrix. The partitions can then be plucked from the main covariance matrix. Here is some code to get $$S$$, $$Sxx$$, $$Syx$$.

> #Read in Data
> names(mydf)<-c("y1", "y2", "x1", "x2")
> #Change column order so S is created in the same order as in OP's question
> mydf2<-data.frame(mydf$$x1, mydf$$x2, mydf$$y1, mydf$$y2)
> names(mydf2)<-c("x1", "x2", "y1", "y2")
> #Print to compare to Table 3.8 in book
x1  x2  y1  y2
1 179 145 191 155
2 201 152 195 149
3 185 149 181 148
4 188 149 183 153
5 171 142 176 144
6 192 152 208 157
> #Obtain full variance covariance matrix
> S <- cov(mydf2)
> #Obtain covariance partitioned matrices
> Sxx <- S[1:2,1:2]
> Syy <- S[3:4,3:4]
> Syx <- Sxy <- S[1:2, 3:4]
> S
x1       x2       y1       y2
x1 100.80667 56.54000 69.66167 51.31167
x2  56.54000 45.02333 46.11167 35.05333
y1  69.66167 46.11167 95.29333 52.86833
y2  51.31167 35.05333 52.86833 54.36000
> Sxx
x1       x2
x1 100.8067 56.54000
x2  56.5400 45.02333
> Syy
y1       y2
y1 95.29333 52.86833
y2 52.86833 54.36000
> Syx
y1       y2
x1 69.66167 51.31167
x2 46.11167 35.05333


We could use your method, which is mathematically correct, to obtain $$Syx$$ too, but it's more work. Here are the alternate calculations:

> #Your method
> n<-25
> xbar<-apply(mydf2, 2, mean)[1:2]
> ybar<-apply(mydf2, 2, mean)[3:4]
> Syx2<-(t(mydf2[1:2])%*%as.matrix(mydf2[3:4])-n*(xbar%*%t(ybar)))/(n-1)
> Syx2
y1       y2
x1 69.66167 51.31167
x2 46.11167 35.05333
>


Note that the $$Syx$$ here matches the computations I provided previously for $$Syx$$.

• Great, thank you! I'm not sure what I did wrong with my implementation. Maybe a sorting issue? Anyway, thank you for your answer and for satisfying my curiosity Feb 25, 2019 at 18:50