Mercurial > repos > nicolas > oghma

comparison computeR2.R @ 28:40664d2d295f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author nicolas
date Tue, 25 Oct 2016 14:38:15 -0400
parents 61780259e2fb
children
comparison
equal deleted inserted replaced
27:203872c77910 28:40664d2d295f
1 ######################################################## 1 ########################################################
2 # 2 #
3 # creation date : 27/06/16 3 # creation date : 27/06/16
4 # last modification : 27/06/16 4 # last modification : 22/10/16
5 # author : Dr Nicolas Beaume 5 # author : Dr Nicolas Beaume
6 # owner : IRRI 6 # owner : IRRI
7 # 7 #
8 ######################################################## 8 ########################################################
9 9
10 log <- file(paste(getwd(), "log_computeR2.txt", sep="/"), open = "wt") 10 # compute r2 by computing the classic formula
11 sink(file = log, type="message") 11 # compare the sum of square difference from target to prediciton
12 12 # to the sum of square difference from target to the mean of the target
13 library("miscTools") 13 computeR2 <- function(target, prediction) {
14 library(randomForest) 14 sst <- sum((target-mean(target))^2)
15 15 ssr <- sum((target-prediction)^2)
16 computeR2 <- function(phenotype, prediction) { 16 return(1-ssr/sst)
17 return(rSquared(phenotype, (phenotype - prediction))[1,1])
18 } 17 }
19 ############################ main ############################# 18 ############################ main #############################
19 # extract argument
20 cmd <- commandArgs(trailingOnly = T) 20 cmd <- commandArgs(trailingOnly = T)
21 source(cmd[1]) 21 source(cmd[1])
22 # load target and prediction
22 phenotype <- read.table(phenotype, sep="\t", h=T)[,1] 23 phenotype <- read.table(phenotype, sep="\t", h=T)[,1]
23 predicted <- read.table(predicted, sep = "\t", h=T)[,2] 24 predicted <- read.table(predicted, sep = "\t", h=T)[,2]
25 # compute r2
24 cat(computeR2(phenotype, predicted), file=out) 26 cat(computeR2(phenotype, predicted), file=out)