prs_nl.r

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),])