4 may 2012

Practicing Script with “ R”: Monitor

These are  samples analyzed by a reference method (column: Protein) and by an analytical method with a certain model (column: IFTpro). The idea is to create a Monitor Report for some basic statistics (RMSEP, Bias, SEP, R,RSQ) to see how well the model performs.
Sample Protein  IFTpro
3      12.85    12.95
4      12.68    12.59
5      11.94    12.12
6      12.07    12.25
7      12.53    12.35
8      11.82    12.20
9      12.58    12.18
10     12.35    12.27
11     12.38    12.32
12     12.15    12.31
13     12.75    12.28
14     12.51    12.07
15     11.92    12.20
16     12.14    12.24
17     12.33    12.27
18     12.15    12.10
20     11.82    11.94
21     11.82    12.05
22     12.36    12.05
23     12.06    11.91
24     11.87    11.98
25     11.81    11.80
26     11.53    11.64
27     11.75    11.84
I take this as a practice with R to write some script.
This is the script:
monitor2<-function(x,y){
 n<-length(y)
 res<-y-x
 par(mfrow=c(2,2))
 hist(res,col="blue")   
 plot(x~y,xlab="predicted",ylab="reference",lty=1)
 abline(0,1,col="blue")
 l<-seq(1:n)
 plot(res~l)
 abline(0,0,col="blue")
 {rmsep<-sqrt(sum((y-x)^2)/n)
 cat("RMSEP:",rmsep,"\n")}
 {(bias<-mean(res))
 cat("Bias :",bias,"\n")}
 {sep<-sd(res)
 cat("SEP  :",sep,"\n")}
 {r<-cor(x,y)
 cat("Corr :",r,"\n")}
 {rsq<-(r^2)
 cat("RSQ  :",rsq,"\n")}   

 }The statistics for this case are:
> monitor2(semola1$Protein,semola1$IFTpro)
RMSEP: 0.2219797
Bias : -0.01083333
SEP  : 0.2264838
Corr : 0.772607
RSQ  : 0.5969215  

And the plots:

I realized that there are a lot of things to improve. to make this script more robust. So I will continue reading tutorials, R help pages, and posts from R blogger,...looking at videos, webinars, reading books,.... to continue improving. Anyway, feel free to take this scrip and add to me feedback.

No hay comentarios:

Publicar un comentario