This is a test version of Biostars. For the public version, visit https://www.biostars.org.
instead of linear regression, do logistic regression plink prs tutorial.

Hi everyone,

Some questions: this is the code on the PLINK tutorial to calculate PRS. https://choishingwan.github.io/PRS-Tutorial/plink/ (go down). This is done in R-studio. The phenotype data I'm working with is not continuous but binary. So 1 stands for control and 2 stands for cases. However here, if I'm understanding correctly, they are working with continuous phenotype. In this script they are doing the linear regression, but since my phenotype is control/case it's binary right. And that means I should be working with logistic regression. However, can someone help me to start?

  library(data.table)
library(magrittr)
p.threshold <- c(0.001,0.05,0.1,0.2,0.3,0.4,0.5)
phenotype <- fread("EUR.height")
pcs <- fread("EUR.eigenvec", header=F) %>%
    setnames(., colnames(.), c("FID", "IID", paste0("PC",1:6)) )
covariate <- fread("EUR.cov")
pheno <- merge(phenotype, covariate) %>%
        merge(., pcs)
null.r2 <- summary(lm(Height~., data=pheno[,-c("FID", "IID")]))$r.squared
prs.result <- NULL
for(i in p.threshold){
    pheno.prs <- paste0("EUR.", i, ".profile") %>%
        fread(.) %>%
        .[,c("FID", "IID", "SCORE")] %>%
        merge(., pheno, by=c("FID", "IID"))

    model <- lm(Height~., data=pheno.prs[,-c("FID","IID")]) %>%
            summary
    model.r2 <- model$r.squared
    prs.r2 <- model.r2-null.r2
    prs.coef <- model$coeff["SCORE",]
    prs.result %<>% rbind(.,
        data.frame(Threshold=i, R2=prs.r2, 
                    P=as.numeric(prs.coef[4]), 
                    BETA=as.numeric(prs.coef[1]),
                    SE=as.numeric(prs.coef[2])))
}
print(prs.result[which.max(prs.result$R2),])
q() # exit R
regression plink prs logistic

1 answer

replace lm with glm(formula, data, family=binomial) (replace formula with your regression formula, and data with your data frame)

Make sure your case control is coded as 0 (control) and 1 (case)

to get the R2, use fmsb::NagelkerkeR2(glm.model)

where glm.model is the glm object.

If these are too complicated, I'd recommend you to use dedicated software such as PRSice-2 to get your PRS results.

Log in to answer this question.