Mercurial > repos > nicolas > oghma

comparison evaluation.R @ 75:dcbee2f7b600 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author nicolas
date Fri, 28 Oct 2016 08:45:47 -0400
parents
children
comparison
equal deleted inserted replaced
74:dea79015dfc9 75:dcbee2f7b600
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)