= R code for Jin's test of a correlation between baseline, X, and change from baseline, Y-X =

(written by Leon Reteig of the Department of Psychology, University of Amsterdam, January 2019)

{{{
corr_change_baseline <- function(X,Y) {
  # Performs a (two-tailed) test for equal variances in the two input vectors. 
  # This is one way to do an unbiased assesment of the correlation between baseline and change from baseline (retest - baseline).
  # In this case, X would be the baseline, and Y the retest.
  #
  # Returns a list of 4 values, following Maloney & Rastogi (1970)
  #   -r: correlation of X-Y vs. X+Y. Assessing significance of this correlation is equivalent to a variance test
  #   -t: t-value constructed from the r-value
  #   -df: degrees of freedom for "r" and "t"
  #   -p: p-value for the variance test. If significant, this means the two vectors do not have equal variances,
  #       suggesting presence of a correlation between baseline and change from baseline.
  #
  # This implementation was inspired by and checked against an Excel spreadsheet, available at:
  # http://imaging.mrc-cbu.cam.ac.uk/statswiki/FAQ/rxxy_correction
  # This page is a great resource for further background on the issue of variance tests / correlations
  #
  # The relevant methods are contained in the following papers:
  #
  #   Jin, P. (1992). Toward a reconceptualization of the law of initial value. Psychological Bulletin 111(1) 176-184.
  #   Maloney, C. J. and Rastogi, S. C. (1970). Significance test for Grubb's estimators. Biometrics 26 671-676.
  #   Myrtek, M. and Foerster, F. (1986). The law of initial value: a rare exception. Biological Psychology 22 227-237
  #   Tu, Y-K. and Gilthorpe, M. S. (2007). Revisiting the relation between change and initial value: A review and evaluation. 
  #     Statistics in Medicine 26 443-457
  
  # This approach of testing the correlation correlation of X-Y and X+Y 
  # is equivalent to testing the variance of X vs. the variance of Y (Tu & Gilthorpe, 2007) 
  r_cor <- cor(X - Y, X + Y)
  n <- length(X) # number of observations 
  df <- n - 2 # degrees of freedom for r and t
  
  t_MR_1970 <- r_cor*sqrt(df) / sqrt(1 - r_cor^2) # t-value proposed by Maloney & Rastogi (1970); equation taken from Tu & Gilthorpe (2007)
  p_MR_1970 <- 2*pt(-abs(t_MR_1970), df) # p-value from t-value (for a two-sided test)
  return(list(r = r_cor, df = df, t = t_MR_1970, p = p_MR_1970))
  
  # Myrtek & Foerster (1986) proposed a different test. The approach is equivalent to Maloney & Rastogi (1970) and yields the same p-value.
  # It is only included here for the sake of completeness.
  
  B_e <- ( var(Y) - var(X) + sqrt((var(Y) - var(X))^2 + 4*(cov(X,Y))^2) ) / (2*cov(X,Y)) # structural relation, from Jin (1992)
  # N.B. this is a simplifaction of the equation in the paper, as "cor(X,Y) * sd(X) * sd(Y)"
  # simply reduces to "cov(X,Y)". The formula as literally in Jin (1992):
  # B_e_paper <- ( var(Y) - var(X) + sqrt((var(Y) - var(X))^2 + 4*(cor(X,Y) * sd(X) * sd(Y))^2) ) / (2*cor(X,Y) * sd(X) * sd(Y))
  # This simplification is also used in the formula below.
  
  # t-value proposed by Myrtek & Foerster, 1986; equation taken from from Jin (1992)
  t_MF_1986 <- sqrt( df * sin(2*(atan(B_e) - atan(1)))^2 *
                       ( (1/4) * (var(X) - var(Y))^2 + cov(X,Y)^2 ) / 
                       ( var(X) * var(Y) - cov(X,Y)^2) )
  p_MF_1986 <- 2*pt(-abs(t_MF_1986), df) #p-value from t-value (for a two-sided test)
  # return(list(df = df, t = t_MF_1986, p = p_MF_1986)) # the t- and p-values are identical to t_MR_1970 and p_MR_1970, so only return those
}
}}}

__Example using the above R function (given in spreadsheet on preceding page)__

{{{
library(foreign)

> dat2 <- data.frame(dat)
> X <- dat2[,1]
> Y <- dat2[,2]
> X
 [1]  1.05842  2.43411 -0.51272  0.96959  4.00000  5.00000  2.00000  3.00000  4.00000  5.00000  6.00000  7.00000  3.74174
[14] -1.12423  1.67558  2.57155  2.81140  4.79701 -0.31238  3.12158  2.45195  2.07938  5.03972  1.72247  1.06231 -0.99752
[27]  2.66932 -4.12602  2.81685  2.97308  3.13575  0.46938  0.46938
> Y
 [1]  3.41853  0.87283  2.76419  5.36829  2.15316  3.47944 -0.66727 -0.07276 -0.06596  1.54776  2.52074 -0.52881  5.48375
[14]  0.59687 -0.56258  2.86513  2.47741 -0.01483  4.12124  1.63788  0.88078  4.90384  3.36825  2.04285  3.68688  5.04538
[27]  1.52683 -1.26871  0.83673  3.58103  1.03052 -1.08757 -1.08757
}}}

This yields the same answer as the spreadsheet

{{{
> corr_change_baseline(X,Y)
$r
[1] 0.1296376

$df
[1] 31

$t
[1] 0.7279344

$p
[1] 0.4721166
}}}