1 Objective

In this practical, we evaluate prediction accuracy for quantitative and binary traits using several statistics. Examples use simple simulated data in R, and we then apply the same code to SBayesRC PGS for the simulated trait with h2 = 0.5 we had analysed in the previous practicals.

Note: R code is shown with a light blue background, while terminal commands are shown with a light orange background.

2 Prediction accuracy (R-square) and bias for quantitative traits

### simulate the data
set.seed(100)  # set this seed so you can get the same result
n = 10000
### quantitative trait
yhat = rnorm(n)
age = as.integer(runif(n, 20, 60))
sex = rbinom(n, 1, 0.55)
y = yhat + scale(age)*sqrt(0.02) +  rnorm(n)*sqrt(0.88)
ids = paste("indi", 1:n, sep="")
data = data.frame(ID=ids, pheno=y, pgs=yhat, age=age, sex=sex)
### linear regression
lmR = lm(pheno ~ age + sex, data=data)  ### reduced module
lmF = lm(pheno ~ age + sex + pgs, data=data) ### full module
### look at the summary results for the two models
summary(lmR)
summary(lmF)
### incremental r-square
summary(lmF)$"r.square" - summary(lmR)$"r.square"
## [1] 0.5196431

The incremental \(R^2\) is often reported as prediction accuracy, quantifying the proportion of phenotypic variance explained by the PGS after adjusting for age and sex.

Prediction bias is assessed from the regression slope for PGS: a slope of 1 indicates no bias.

3 Nagelkerke’s R-square for diseases (binary traits)

### simulate the data
### binary trait
set.seed(101)
yhat = rnorm(n)
age = as.integer(runif(n, 20, 60))
sex = rbinom(n, 1, 0.55)
y = yhat + scale(age)*sqrt(0.02) +  rnorm(n)*sqrt(0.88)
y = as.numeric(y>=0)
ids = paste("indi", 1:n, sep="")
data = data.frame(ID=ids, pheno=y, pgs=yhat, age=age, sex=sex)
### logistic regression
glmR = glm(pheno ~ age + sex, data=data, family=binomial(logit))  ### reduced module
glmF = glm(pheno ~ age + sex + pgs, data=data, family=binomial(logit))  ### full module
### look at the summary results for the two models
summary(glmR)
summary(glmF)
### log-likelihood
N = nrow(data)
LLF = logLik(glmF)
LLR = logLik(glmR)

### Cox&Snell R2
CSv <-  1-exp((2/N)*(LLR[1]-LLF[1]))
CSv
## [1] 0.3484616
### Nagelkerke's R2
NKv <- CSv/(1-exp((2/N)*LLR[1]))
NKv
## [1] 0.4653179

Nagelkerke’s \(R^2\) is a pseudo-\(R^2\) for logistic regression that approximates the proportion of variance explained on the observed binary scale.

4 AUC

### AUC
### install.packages('pROC')
library('pROC')
aucF = auc(data$pheno, glmF$linear.predictors)  ### AUC for full module
aucR = auc(data$pheno, glmR$linear.predictors)  ### AUC for reduced module
aucF; aucR
## Area under the curve: 0.8503
## Area under the curve: 0.5387
aucF - aucR
## [1] 0.311632
### draw the ROC
#install.packages("PredictABEL")
library(PredictABEL)
plotROC(data=data, cOutcome=2, predrisk=glmF$linear.predictors)

## AUC [95% CI] for the model 1 :  0.85 [ 0.843  -  0.858 ]

AUC measures the discriminative ability of the model. It quantifies the probablity of ranking a case higher than a control based on the predictor. The difference aucF - aucR shows how much PGS improves case/control classification beyond age and sex.

5 Variance explained on liability scale

### function to convert R-square from 0-1 observed scale
### to liability scale
h2l_R2 <- function(k, r2, p) {
  # k baseline disease risk
  # r2 from a linear regression model of genomic profile risk score
  # p proportion of sample that are cases
  # calculates proportion of variance explained on the liability scale
  # from ABC at http://www.complextraitgenomics.com/software/
  # Lee SH, Goddard ME, Wray NR, Visscher PM. (2012) A better coefficient of determination for genetic profile analysis. Genet Epidemiol. 2012 Apr;36(3):214-24.
  x = qnorm(1-k)
  z = dnorm(x)
  i = z/k
  C = k*(1-k)*k*(1-k)/(z^2*p*(1-p))
  theta = i*((p-k)/(1-k))*(i*((p-k)/(1-k))-x)
  h2l_R2 = C*r2 / (1 + C*theta*r2)
  return(h2l_R2)
}
K=0.1  ## population prevalence
P = sum(data$pheno)/nrow(data)  ### proportion of cases in the sample
P
## [1] 0.5029
### linear regression with 0/1 values
lmR = lm(pheno~age+sex, data=data)
lmF = lm(pheno~age+sex+pgs, data=data)
#R2v = summary(lmF)$"r.square" - summary(lmR)$"r.square"
#To strictly follow the equation 7 in Lee et al.
R2v = 1-exp((2/N)*(logLik(lmR)[1]-logLik(lmF)[1]))
R2v
## [1] 0.3459491
### convert to liability scale
h2l_R2(K,R2v,P)
## [1] 0.4240657

Converts observed-scale \(R^2\) to liability-scale, accounting for case-control ascertainment and population prevalence. The liability-scale \(R^2\) expresses the variance explained on an underlying continuous liability (risk) for disease, adjusting for population prevalence so it is comparable to heritability on the liability scale. For many diseases the liability is a latent (unobserved) variable, but for some diseases the liability can be measured directly, for example, body mass index (BMI) can be interpreted as the liability for obesity.

6 Decile Odds Ratio

### cut into deciles
data$decile = cut(data$pgs, breaks=c(quantile(data$pgs, probs=seq(0,1,by=0.1))), labels=1:10, include.lowest=T)
### calculate manually the odds in each decile
#### install.packages("tidyverse")
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
data %>% group_by(decile) %>% 
    summarise(n_case=sum(pheno==1), n_control=sum(pheno==0)) %>%
    mutate(odds = n_case/n_control) %>%
    mutate(ORs = odds/odds[1])
## # A tibble: 10 × 5
##    decile n_case n_control    odds    ORs
##    <fct>   <int>     <int>   <dbl>  <dbl>
##  1 1          46       954  0.0482   1   
##  2 2         133       867  0.153    3.18
##  3 3         241       759  0.318    6.59
##  4 4         332       668  0.497   10.3 
##  5 5         462       538  0.859   17.8 
##  6 6         546       454  1.20    24.9 
##  7 7         668       332  2.01    41.7 
##  8 8         759       241  3.15    65.3 
##  9 9         870       130  6.69   139.  
## 10 10        972        28 34.7    720.
### calculate ORs using logistic regression
glmD <- glm(pheno ~ decile, data = data, family = binomial(logit))
### Odds for being a case compared to control in each decile
ORD <- exp(glmD$coefficients)
ORD
##  (Intercept)      decile2      decile3      decile4      decile5      decile6 
##   0.04821803   3.18143523   6.58515209  10.30747201  17.80943915  24.94177361 
##      decile7      decile8      decile9     decile10 
##  41.72812991  65.31535270 138.79264214 719.94409936
### Plot odds
ORDL <- exp(glmD$coefficients-1.96*summary(glmD)$coefficients[,2])
ORDH <- exp(glmD$coefficients+1.96*summary(glmD)$coefficients[,2])
plot(ORD,ylim=c(min(ORDL),max(ORDH)))
arrows(seq(1,10,1), ORD, seq(1,10,1), ORDH, angle=90,length=0.10) # Draw error bars
## Warning in arrows(seq(1, 10, 1), ORD, seq(1, 10, 1), ORDH, angle = 90, length =
## 0.1): zero-length arrow is of indeterminate angle and so skipped
arrows(seq(1,10,1), ORD, seq(1,10,1), ORDL, angle=90,length=0.10) # Draw error bars
## Warning in arrows(seq(1, 10, 1), ORD, seq(1, 10, 1), ORDL, angle = 90, length =
## 0.1): zero-length arrow is of indeterminate angle and so skipped

This plot shows how disease odds change across higher PGS deciles. ORs relative to the lowest decile illustrate risk stratification ability.

7 Using PGS from SBayesRC

We have generated PGS from SBayesRC for the simulated trait with ~300K SNPs in the target sample. Recall we have

  • Phenotype file: target_phenotypes.txt
  • Covariate file: target_covariates.txt
  • PGS file: target_sbayesrc.profile

Read these files into R.

library(data.table)
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## The following object is masked from 'package:base':
## 
##     %notin%
qt <- fread("target_phenotypes.txt") %>%
  setnames(., names(.), c("FID", "IID", "pheno"))
covar <- fread("target_covariates.txt")
prs <- fread("target_sbayesrc.profile") %>%
  select(FID, IID, SCORESUM) %>%
  setnames(., names(.)[3], 'SBRC_SCORE')
qt <- merge(qt, prs, by=c("FID", "IID")) %>%
  merge(., covar, by=c("FID", "IID"))

Based on liability theory, we convert quantitative phenotypes to binary disease status assuming population prevalence K. The code below performs this conversion with K = 0.1.

K <- 0.1 # disease prevalence
t <- qnorm(1-K)

bt <- qt %>%
  mutate(pheno_scaled = scale(pheno), 
         disease_status = as.numeric(pheno_scaled > t)) %>%
  select(FID, IID, disease_status)

table(bt$disease_status)
## 
##   0   1 
## 442  52

8 Distribution of PRS

library(ggplot2)
# add PRS and covariates
bt <- merge(bt, prs, by=c("FID", "IID")) %>%
  merge(., covar, by=c("FID", "IID"))

plt <- ggplot(data = bt, aes(x = SBRC_SCORE, fill = as.factor(disease_status))) +
  geom_density(alpha=0.5) +
  scale_fill_discrete(labels = c("controls", "cases")) +
  labs(title = "Distribution of PRS", 
       x = "PRS - SBayesRC", 
       y = "Density", 
       fill = "Disease status") +
  theme(plot.title = element_text(hjust=0.5))
plt

Q1: Why do PRS of cases and controls overlap as shown in the graph?

Major considerations
  • PRS captures only part of the genetic risk
  • Environmental and lifestyle factors matter
  • PRS itself is an estimate with noise depending on GWAS power

Now try to compute the statistics above in the target sample. Write your own code before checking the provided solution.

Q2: What’s the prediction accuracy and bias for quantitative phenotypes?

Show code 
# linear regression
covar_cols <- setdiff(names(covar), c("FID", "IID"))
null_covar <- reformulate(covar_cols, response = "pheno")
full_sbrc <- reformulate(c(covar_cols, "SBRC_SCORE"), response = "pheno")

lmR <- lm(null_covar, data = qt) # reduced model without PRS
lmF <- lm(full_sbrc, data = qt) # full model

# check results of the two models
summary(lmR)
## 
## Call:
## lm(formula = null_covar, data = qt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -45.670  -8.182   0.334   8.977  39.732 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.04219    1.92678   0.022   0.9825  
## SEX           0.40850    1.20819   0.338   0.7354  
## PC1         -25.70002   13.32124  -1.929   0.0543 .
## PC2          13.67180   13.27758   1.030   0.3037  
## PC3          -8.82483   13.27675  -0.665   0.5066  
## PC4         -11.22673   13.27590  -0.846   0.3982  
## PC5           9.86101   13.27984   0.743   0.4581  
## PC6          12.83185   13.28308   0.966   0.3345  
## PC7          -2.48914   13.27798  -0.187   0.8514  
## PC8           3.03224   13.31541   0.228   0.8200  
## PC9          -0.18922   13.29505  -0.014   0.9887  
## PC10         -5.95964   13.29891  -0.448   0.6543  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.28 on 482 degrees of freedom
## Multiple R-squared:  0.01583,    Adjusted R-squared:  -0.006626 
## F-statistic: 0.705 on 11 and 482 DF,  p-value: 0.7342
summary(lmF)
## 
## Call:
## lm(formula = full_sbrc, data = qt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.290  -8.037   0.184   7.978  34.779 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.28341    1.68632  -1.947   0.0521 .  
## SEX           0.64226    1.04493   0.615   0.5391    
## PC1         -16.60865   11.54126  -1.439   0.1508    
## PC2           9.38372   11.48649   0.817   0.4144    
## PC3          -1.01487   11.49711  -0.088   0.9297    
## PC4          -9.68393   11.48078  -0.843   0.3994    
## PC5          12.82078   11.48588   1.116   0.2649    
## PC6           8.98699   11.49028   0.782   0.4345    
## PC7          -4.75802   11.48331  -0.414   0.6788    
## PC8           4.72167   11.51507   0.410   0.6820    
## PC9           4.25937   11.50197   0.370   0.7113    
## PC10          0.80858   11.51221   0.070   0.9440    
## SBRC_SCORE    0.71555    0.05595  12.790   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.48 on 481 degrees of freedom
## Multiple R-squared:  0.2656, Adjusted R-squared:  0.2473 
## F-statistic:  14.5 on 12 and 481 DF,  p-value: < 2.2e-16


Q3: What’s the Nagelkerke’s R-square for the binary phenotypes?

Show code 
# logistic regression
null_covar <- reformulate(covar_cols, response = "disease_status")
full_sbrc <- reformulate(c(covar_cols, "SBRC_SCORE"), response = "disease_status")

glmR <- glm(null_covar, data = bt, family=binomial(logit))
glmF <- glm(full_sbrc, data = bt, family=binomial(logit))
summary(glmR)
## 
## Call:
## glm(formula = null_covar, family = binomial(logit), data = bt)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.37277    0.49412  -4.802 1.57e-06 ***
## SEX           0.06263    0.30023   0.209    0.835    
## PC1          -0.93421    3.22774  -0.289    0.772    
## PC2          -6.98759    4.55295  -1.535    0.125    
## PC3           2.79053    3.26495   0.855    0.393    
## PC4         -13.53058    8.33185  -1.624    0.104    
## PC5          -0.87038    5.75304  -0.151    0.880    
## PC6           3.34609    4.48202   0.747    0.455    
## PC7           5.23742    4.73470   1.106    0.269    
## PC8          -5.55108    4.28843  -1.294    0.196    
## PC9           1.34317    3.77367   0.356    0.722    
## PC10         -3.06575    3.63055  -0.844    0.398    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 332.46  on 493  degrees of freedom
## Residual deviance: 324.23  on 482  degrees of freedom
## AIC: 348.23
## 
## Number of Fisher Scoring iterations: 6
summary(glmF)
## 
## Call:
## glm(formula = full_sbrc, family = binomial(logit), data = bt)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.43658    0.58440  -5.880 4.09e-09 ***
## SEX           0.08110    0.33012   0.246    0.806    
## PC1           0.67681    3.53695   0.191    0.848    
## PC2          -6.97004    4.89502  -1.424    0.154    
## PC3           4.30244    3.59584   1.197    0.231    
## PC4         -14.25321    9.08681  -1.569    0.117    
## PC5          -1.76492    6.01188  -0.294    0.769    
## PC6           4.31680    4.85589   0.889    0.374    
## PC7           3.08873    4.69878   0.657    0.511    
## PC8          -3.95919    4.69095  -0.844    0.399    
## PC9           3.61984    4.29468   0.843    0.399    
## PC10         -0.97340    3.84437  -0.253    0.800    
## SBRC_SCORE    0.12924    0.02002   6.457 1.07e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 332.46  on 493  degrees of freedom
## Residual deviance: 271.63  on 481  degrees of freedom
## AIC: 297.63
## 
## Number of Fisher Scoring iterations: 7
# log-likelihood
N= nrow(bt)
LLF = logLik(glmF)
LLR = logLik(glmR)
c(LLF, LLR)
## [1] -135.8130 -162.1168
# Nagelkerke's R2
NKv <- CSv / (1-exp((2/N)*LLR[1]))
NKv
## [1] 0.724068


Q4: What’s the AUC for the binary phenotypes?

Show code 
# add PRS and covariates
aucF <- auc(bt$disease_status, glmF$linear.predictors)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
aucR <- auc(bt$disease_status, glmR$linear.predictors)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
c(aucF, aucR, aucF-aucR)
## [1] 0.8047337 0.6493648 0.1553690
plotROC(data = bt, cOutcome = 3, predrisk=glmF$linear.predictors)

## AUC [95% CI] for the model 1 :  0.805 [ 0.747  -  0.863 ]


Q5: What’s the R-square on liability scale for the binary phenotypes?

Show code 
K=0.1  ## population prevalence
P = sum(bt$disease_status)/nrow(bt)  ### proportion of cases in the sample
P
## [1] 0.1052632
### linear regression with 0/1 values
lmR = lm(null_covar, data=bt)
lmF = lm(full_sbrc, data=bt)
#R2v = summary(lmF)$"r.square" - summary(lmR)$"r.square"
#To strictly follow the equation 7 in Lee et al.
R2v = 1-exp((2/N)*(logLik(lmR)[1]-logLik(lmF)[1]))
R2v
## [1] 0.1025564
### convert to liability scale
h2l_R2(K,R2v,P)
## [1] 0.2874458


Q6: What are the decile odds ratios for the binary phenotypes?

Show code 
bt$decile = cut(bt$SBRC_SCORE, breaks=c(quantile(bt$SBRC_SCORE, probs=seq(0,1,by=0.25))), labels=1:4, include.lowest=T)

### calculate manually the odds in each decile
bt %>% group_by(decile) %>% 
    summarise(n_case=sum(bt==1), n_control=sum(bt==0)) %>%
    mutate(odds = n_case/n_control) %>%
    mutate(ORs = odds/odds[1])
## # A tibble: 4 × 5
##   decile n_case n_control  odds   ORs
##   <fct>   <int>     <int> <dbl> <dbl>
## 1 1         415       442 0.939     1
## 2 2         415       442 0.939     1
## 3 3         415       442 0.939     1
## 4 4         415       442 0.939     1
### calculate ORs using logistic regression
glmD <- glm(disease_status ~ decile, data = bt, family = binomial(logit))
### Odds for being a case compared to control in each decile
ORD <- exp(glmD$coefficients)
ORD
##  (Intercept)      decile2      decile3      decile4 
## 3.181005e-09 2.481837e+07 3.398552e+07 1.047887e+08
### Plot odds
ORDL <- exp(glmD$coefficients - 1.96 * summary(glmD)$coefficients[,2])
ORDH <- exp(glmD$coefficients + 1.96 * summary(glmD)$coefficients[,2])
# Ensure finite limits and use the correct index length for plotting
# Plot points and draw error bars only where CIs are finite
vals <- c(ORD, ORDL, ORDH)
if (!any(is.finite(vals))) {
  cat('Cannot plot ORs: all values are non-finite\n')
} else {
  ylim <- range(vals[is.finite(vals)])
  plot(ORD, ylim = ylim)
  idx <- seq_along(ORD)
  finiteCI <- is.finite(ORDL) & is.finite(ORDH)
  if (any(finiteCI)) {
    arrows(idx[finiteCI], ORD[finiteCI], idx[finiteCI], ORDH[finiteCI], angle = 90, length = 0.10)
    arrows(idx[finiteCI], ORD[finiteCI], idx[finiteCI], ORDL[finiteCI], angle = 90, length = 0.10)
  }
}

Note that here we use the median decile as the reference as there is zero case in the first decile.