Created
June 7, 2021 16:49
-
-
Save nievergeltlab/d2565c5135099da84c5a1b4851b905a3 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(MASS) | |
library(data.table) | |
require(foreign) | |
require(nnet) | |
require(ggplot2) | |
require(reshape2) | |
library(lmtest) | |
#if it complains about this package: | |
#install.packages('zoo',repos="https://cran.r-project.org") | |
#install.packages('lmtest',repos="https://cran.r-project.org") | |
#Load pcs | |
pcs <- read.table('phenotypes_covariates/mrs2_eur_mrs_cov_vizeli.txt',header=T,na.strings=c("NA","-9")) # | |
#Load phenotypes and covariates | |
pheno <- read.table('phenotypes_covariates/mrs2_eur_phenos_vizeli.txt',header=T,na.strings=c("NA","-9")) | |
#Load PRS | |
prs <- read.table('prs_outputs/ptsd_prs/freeze2_prs_acq.LastHalfAcqDiscrim.best',header=T,na.strings=c("NA","-9")) | |
#index of PRS names, to loop over all prs | |
prs_scores <- names(prs)[-c(1,2)] | |
#Notice that all.x=T, because the family file length must be preserved! | |
d1 <- merge(pcs,pheno,by=c("FID","IID"),all.x=T,suffixes=c("_fam","")) | |
dm <- merge(d1,prs,by=c("FID","IID"),all.x=T,suffixes=c("_cov","")) | |
dm <- subset(dm,!is.na(dm[,prs_scores[1]])) #just filter out anyone without a score, this will make model likelihoods easier.. | |
##Model for latent class | |
#Make a matrix for outputs | |
outmat <- as.data.frame(matrix(nrow=length(prs_scores),ncol=5)) | |
row.names(outmat) <- prs_scores | |
#null model | |
m0 <- multinom(latent_class ~ C1 + C2 +C3 +C4+C5,data=dm) | |
#loop over all PRS | |
for (prses in prs_scores) | |
{ | |
m1 <- multinom(latent_class ~ dm[,prses] + C1 + C2 +C3 +C4+C5,data=dm) | |
outmat[prses,] <- as.numeric(lrtest(m0,m1)[2,]) | |
} | |
names(outmat) <-c("DF","LogLik","NPARM","Chisq","P") | |
write.table(outmat,file="PRS_latent_class.txt") | |
#Then just sort this file to find the best p-value threshold | |
head(outmat[order(outmat$P),]) | |
### for suvival analysis | |
library(survival) | |
dm$Surv_ex <- Surv(time=dm$extinction_time,event=dm$survival_status) | |
#Make a matrix for outputs | |
outmat <- as.data.frame(matrix(nrow=length(prs_scores),ncol=5)) | |
row.names(outmat) <- prs_scores | |
#null model | |
m0 <- coxph(Surv_ex ~ C1 + C2 +C3 +C4+C5,data=dm) | |
#loop over all PRS | |
for (prses in prs_scores) | |
{ | |
m1 <- coxph(Surv_ex ~ dm[,prses] + C1 + C2 +C3 +C4+C5,data=dm) | |
outmat[prses,] <- as.numeric(lrtest(m0,m1)[2,]) | |
} | |
names(outmat) <-c("DF","LogLik","NPARM","Chisq","P") | |
write.table(outmat,file="PRS_surv.txt") | |
#Then just sort this file to find the best p-value threshold | |
head(outmat[order(outmat$P),]) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment