This analysis was to compare different methods for estimating the statistical significance of prediction differences across polygenic scoring methods, and models in general.
Here I simulate data with one continuous outcome, and two correlated continuous sets of predictions The correlation between the predictors and the outcome is varied, and the correlation between the predictors themselves is varied.
As I do in my other analyses, I estimate the correlation between the predictors and the outcome to determine their predictive utility. Then I compare different methods by comparing the observed-predicted correlations. Several method comparing correlation are considered:
- Two-sample Z-test
- Compares estimates from two populations and does not account for the correlation between predictors
- Permutation based
- Randomises the phenotype, retaining the correlation between predictors, and then estimates the number times the difference in correlation is larger than the observed difference.
- Cox test
- Method for comparing non-nested models
- Pearson and Filon’s (1898):
- Method for comparing correlations between variables that are correlated and measured in a single sample
- Implemented using cocor.dep.groups.overlap function in the cocor package
- Other methods implemented by this function are highly concordant
- Fisher r-to-z :
- Method for comparing correlations between variables that are independent and measures in different samples.
- Implemented using psych package
- Accounts for non-normal error of correlations
- Williams test :
- Method for comparing correlations between variables that are dependent and measured in the same sample.
- Implemented using psych package
- Accounts for non-normal error of correlations
We will simulate the following scenarios:
- Estimates from two independent samples
- Estimates from one sample, but predictors are uncorrelated
- Estimates from one sample, and predictors are highly correlated
Estimates from two independent samples
library(data.table)
set.seed(1)
# Sample 1
N<-200
y1<-rnorm(N)
x1<-scale(y1+rnorm(N,0,3))
dat1<-data.table(y1,x1)
cor(dat1)
## y1 V1
## y1 1.0000000 0.2743249
## V1 0.2743249 1.0000000
# Sample 2
N<-200
y2<-rnorm(N)
x2<-scale(y2+rnorm(N,0,10))
dat2<-data.table(y2,x2)
cor(dat2)
## y2 V1
## y2 1.00000000 0.06274804
## V1 0.06274804 1.00000000
### Derive models using each predictor
mod1<-lm(y1 ~ x1, data=dat1)
mod2<-lm(y2 ~ x2, data=dat2)
dat1$pred1<-predict(mod1, dat1)
dat2$pred2<-predict(mod2, dat2)
dat_cor<-data.frame(model1='x1',
model2='x2',
cor_x1_x2=NA,
cor_y_x1=coef(summary(mod1))[2,1],
cor_y_x1_se=coef(summary(mod1))[2,2],
cor_y_x2=coef(summary(mod2))[2,1],
cor_y_x2_se=coef(summary(mod2))[2,2],
cor_diff=coef(summary(mod1))[2,1]-coef(summary(mod2))[2,1])
### Test difference between predictors using a Z-test
dat_cor$cor_diff_Ztest_P<-pnorm(-(dat_cor$cor_diff/sqrt((dat_cor$cor_y_x1_se^2)+(dat_cor$cor_y_x2_se^2))))
### Test difference between predictors using a permutation test
n_perm<-500
diff<-NULL
for(i in 1:n_perm){
y1_sample<-sample(y1)
y2_sample<-sample(y2)
mod1_i<-lm(y1_sample ~ x1, data=dat1)
mod2_i<-lm(y2_sample ~ x2, data=dat2)
diff_i<-coef(summary(mod1_i))[2,1]-coef(summary(mod2_i))[2,1]
diff<-c(diff,diff_i)
}
dat_cor$cor_diff_perm_P<-sum(diff > dat_cor$cor_diff[1])/n_perm
### Test difference between predictors using coxtest
# Not possible as different samples
### Test difference between predictors using Fisher's Z transformation
library(cocor)
dat_cor$cocor_fisherZ_P<-cocor.indep.groups(r1.jk=dat_cor$cor_y_x1, r2.hm=dat_cor$cor_y_x2, n1=N, n2=N, alternative='greater')@fisher1925$p.value
library(psych)
dat_cor$psych_fisherZ_P<-paired.r(xy=dat_cor$cor_y_x1, xz=dat_cor$cor_y_x2, n=N, n2=N, twotailed=F)$p
dat_cor
## model1 model2 cor_x1_x2 cor_y_x1 cor_y_x1_se cor_y_x2 cor_y_x2_se
## 1 x1 x2 NA 0.2548745 0.06349504 0.06714255 0.07589417
## cor_diff cor_diff_Ztest_P cor_diff_perm_P cocor_fisherZ_P psych_fisherZ_P
## 1 0.1877319 0.02890093 0.046 0.02747978 0.02747978
The Fisher’s r-to-z transformation is a commonly used approach for comparing correlations from different samples. The two-sample z-test does not account for the non-normal error distribution of Pearson correlations. Results indicate the fisher r-to-z method is concordant with the two-sample z-test. The permutation-based approach is more conservative than other methods.
Estimates from one sample, but predictors are uncorrelated
library(data.table)
set.seed(1)
N<-200
y<-rnorm(N)
x1<-as.numeric(scale(y+rnorm(N,0,3)))
x2<-as.numeric(scale(y+rnorm(N,0,10)))
dat<-data.table(y,x1,x2)
cor(dat)
## y x1 x2
## y 1.0000000 0.27432489 0.15351734
## x1 0.2743249 1.00000000 0.01848117
## x2 0.1535173 0.01848117 1.00000000
### Derive models using each predictor
mod1<-lm(y ~ x1, data=dat)
mod2<-lm(y ~ x2, data=dat)
dat$pred1<-predict(mod1, dat)
dat$pred2<-predict(mod2, dat)
dat_cor<-data.frame(model1='x1',
model2='x2',
cor_x1_x2=cor(dat$x1,dat$x2),
cor_y_x1=coef(summary(mod1))[2,1],
cor_y_x1_se=coef(summary(mod1))[2,2],
cor_y_x2=coef(summary(mod2))[2,1],
cor_y_x2_se=coef(summary(mod2))[2,2],
cor_diff=coef(summary(mod1))[2,1]-coef(summary(mod2))[2,1])
### Test difference between predictors using a Z-test
dat_cor$cor_diff_Ztest_P<-pnorm(-(dat_cor$cor_diff/sqrt((dat_cor$cor_y_x1_se^2)+(dat_cor$cor_y_x2_se^2))))
### Test difference between predictors using a permutation test
n_perm<-500
diff<-NULL
for(i in 1:n_perm){
y_sample<-sample(y)
mod1_i<-lm(y_sample ~ x1, data=dat)
mod2_i<-lm(y_sample ~ x2, data=dat)
diff_i<-coef(summary(mod1_i))[2,1]-coef(summary(mod2_i))[2,1]
diff<-c(diff,diff_i)
}
dat_cor$cor_diff_perm_P<-sum(diff > dat_cor$cor_diff[1])/n_perm
### Test difference between predictors using coxtest
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
dat_cor$cox_diff_P<-coxtest(mod1, mod2)$P[2]
### Test difference between predictors using Pearson
library(cocor)
dat_cor$pearson_diff_P<-cocor.dep.groups.overlap(dat_cor$cor_y_x1, dat_cor$cor_y_x2, dat_cor$cor_x1_x2, N,alternative='greater')@pearson1898$p.value
### Test difference between predictors using Williams's Test
library(psych)
dat_cor$williams_diff_P<-paired.r(xy=dat_cor$cor_y_x1, xz=dat_cor$cor_y_x2, yz=dat_cor$cor_x1_x2, n=N, twotailed=F)$p
dat_cor
## model1 model2 cor_x1_x2 cor_y_x1 cor_y_x1_se cor_y_x2 cor_y_x2_se cor_diff
## 1 x1 x2 0.01848117 0.2548745 0.06349504 0.1426325 0.06524537 0.112242
## cor_diff_Ztest_P cor_diff_perm_P cox_diff_P pearson_diff_P williams_diff_P
## 1 0.1088132 0.112 0 0.1205402 0.1230892
Apart fromt the coxtest method, results are similar across methods, with pearson and williams methods being highly concordant.
Estimates from one sample, and predictors are highly correlated
library(data.table)
set.seed(1)
N<-200
y<-rnorm(N)
set.seed(2)
x1<-as.numeric(scale(y+rnorm(N,0,3)))
set.seed(3)
x2<-as.numeric(scale(x1+rnorm(N,0,1)))
dat<-data.table(y,x1,x2)
cor(dat)
## y x1 x2
## y 1.0000000 0.2813946 0.1928833
## x1 0.2813946 1.0000000 0.6928299
## x2 0.1928833 0.6928299 1.0000000
### Derive models using each predictor
mod1<-lm(y ~ x1, data=dat)
mod2<-lm(y ~ x2, data=dat)
dat$pred1<-predict(mod1, dat)
dat$pred2<-predict(mod2, dat)
dat_cor<-data.frame(model1='x1',
model2='x2',
cor_x1_x2=cor(dat$x1,dat$x2),
cor_y_x1=coef(summary(mod1))[2,1],
cor_y_x1_se=coef(summary(mod1))[2,2],
cor_y_x2=coef(summary(mod2))[2,1],
cor_y_x2_se=coef(summary(mod2))[2,2],
cor_diff=coef(summary(mod1))[2,1]-coef(summary(mod2))[2,1])
### Test difference between predictors using a Z-test
dat_cor$cor_diff_Ztest_P<-pnorm(-(dat_cor$cor_diff/sqrt((dat_cor$cor_y_x1_se^2)+(dat_cor$cor_y_x2_se^2))))
### Test difference between predictors using a permutation test
n_perm<-500
diff<-NULL
for(i in 1:n_perm){
y_sample<-sample(y)
mod1_i<-lm(y_sample ~ x1, data=dat)
mod2_i<-lm(y_sample ~ x2, data=dat)
diff_i<-coef(summary(mod1_i))[2,1]-coef(summary(mod2_i))[2,1]
diff<-c(diff,diff_i)
}
dat_cor$cor_diff_perm_P<-sum(diff > dat_cor$cor_diff[1])/n_perm
### Test difference between predictors using coxtest
library(lmtest)
dat_cor$cox_diff_P<-coxtest(mod1, mod2)$P[2]
### Test difference between predictors using Pearson
library(cocor)
dat_cor$pearson_diff_P<-cocor.dep.groups.overlap(dat_cor$cor_y_x1, dat_cor$cor_y_x2, dat_cor$cor_x1_x2, N,alternative='greater')@pearson1898$p.value
### Test difference between predictors using Williams's Test
library(psych)
dat_cor$williams_diff_P<-paired.r(xy=dat_cor$cor_y_x1, xz=dat_cor$cor_y_x2, yz=dat_cor$cor_x1_x2, n=N, twotailed=F)$p
dat_cor
## model1 model2 cor_x1_x2 cor_y_x1 cor_y_x1_se cor_y_x2 cor_y_x2_se
## 1 x1 x2 0.6928299 0.2614429 0.06336002 0.1792073 0.06478817
## cor_diff cor_diff_Ztest_P cor_diff_perm_P cox_diff_P pearson_diff_P
## 1 0.08223559 0.1820774 0.068 2.185878e-06 0.06303985
## williams_diff_P
## 1 0.06448287
Again the coxtest method is very different from other methods. The two-sample Z test is now deviating from the other methods because it doesn’t account for the correlation between predictors. The results for the permutation, pearson and williams method are highly concordant.
The psych package is a solid package and the Williams method is recommended by Steiger who worked alot in this area. It is also faster than the permutation based approach. From here on use the Williams test to test for significant differences between correlations between outcomes and correlated predictors.