|
75
|
1
|
|
|
2 ########################################################
|
|
|
3 #
|
|
|
4 # creation date : 08/01/16
|
|
|
5 # last modification : 23/06/16
|
|
|
6 # author : Dr Nicolas Beaume
|
|
|
7 # owner : IRRI
|
|
|
8 #
|
|
|
9 ########################################################
|
|
|
10
|
|
|
11 # compute correlation
|
|
|
12 predictionCorrelation <- function(prediction, obs) {
|
|
|
13 return(cor(prediction, obs, method = "pearson"))
|
|
|
14 }
|
|
|
15
|
|
|
16 # compute r2 by computing the classic formula
|
|
|
17 # compare the sum of square difference from target to prediciton
|
|
|
18 # to the sum of square difference from target to the mean of the target
|
|
|
19 r2 <- function(target, prediction) {
|
|
|
20 sst <- sum((target-mean(target))^2)
|
|
|
21 ssr <- sum((target-prediction)^2)
|
|
|
22 return(1-ssr/sst)
|
|
|
23 }
|
|
|
24
|
|
|
25 ################################## main function ###########################
|
|
|
26
|
|
|
27 evaluatePredictions <- function(data, prediction=NULL, traitValue=NULL, doR2=F,
|
|
|
28 folds=NULL, classes=NULL, doCor=F) {
|
|
|
29 eval <- NULL
|
|
|
30 # case of r2
|
|
|
31 if(doR2) {
|
|
|
32 eval <- c(eval, list(r2=NULL))
|
|
|
33 }
|
|
|
34 if(doCor) {
|
|
|
35 eval <- c(eval, list(cor=NULL))
|
|
|
36 }
|
|
|
37 for (i in 1:length(folds)) {
|
|
|
38 train <- unlist(folds[-i])
|
|
|
39 test <- folds[[i]]
|
|
|
40 if(doR2) {
|
|
|
41 if(!is.null(traitValue)) {
|
|
|
42 eval$r2 <- c(eval$r2, r2(traitValue[test], prediction[[i]]))
|
|
|
43 } else {
|
|
|
44 eval$r2 <- c(eval$r2, NA)
|
|
|
45 }
|
|
|
46 }
|
|
|
47 # case of correlation
|
|
|
48 if(doCor) {
|
|
|
49 if(!is.null(traitValue)) {
|
|
|
50 eval$cor <- c(eval$cor, cor(traitValue[test], prediction[[i]]))
|
|
|
51 } else {
|
|
|
52 eval$cor <- c(eval$cor, NA)
|
|
|
53 }
|
|
|
54 }
|
|
|
55 }
|
|
|
56 # print output
|
|
|
57 print(eval)
|
|
|
58 if(doR2) {
|
|
|
59 print(paste("mean r2 : ", mean(eval$r2), "standard deviation :", sd(eval$r2)))
|
|
|
60 }
|
|
|
61 if(doCor) {
|
|
|
62 print(paste("mean correlation : ", mean(eval$cor), "standard deviation :", sd(eval$cor)))
|
|
|
63 }
|
|
|
64 }
|
|
|
65 ############################ main #############################
|
|
|
66
|
|
|
67 # load arguments
|
|
|
68 cmd <- commandArgs(trailingOnly = T)
|
|
|
69 source(cmd[1])
|
|
|
70 # set which evaluation are used
|
|
|
71 if(as.integer(doR2) == 1) {
|
|
|
72 doR2 <- T
|
|
|
73 }
|
|
|
74 if(as.integer(doCor) == 1) {
|
|
|
75 doCor <- T
|
|
|
76 }
|
|
|
77 # load genotype & phenotype
|
|
|
78 con = file(genotype)
|
|
|
79 genotype <- readLines(con = con, n = 1, ok=T)
|
|
|
80 close(con)
|
|
|
81 genotype <- read.table(genotype, sep="\t",h=T)
|
|
|
82 phenotype <- read.table(phenotype, sep="\t",h=T)[,1]
|
|
|
83 # load prediction
|
|
|
84 con = file(prediction)
|
|
|
85 prediction <- readLines(con = con, n = 1, ok=T)
|
|
|
86 close(con)
|
|
|
87 prediction <- readRDS(prediction)
|
|
|
88 # load folds
|
|
|
89 con = file(folds)
|
|
|
90 folds <- readLines(con = con, n = 1, ok=T)
|
|
|
91 close(con)
|
|
|
92 folds <- readRDS(folds)
|
|
|
93 # evaluation
|
|
|
94 evaluatePredictions(data=genotype, prediction=prediction, traitValue=phenotype, doR2=doR2,
|
|
|
95 folds=folds, doCor=doCor) |
