-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplumb_proast.R
153 lines (128 loc) · 5.44 KB
/
plumb_proast.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
# script name
# plumb_proast.R
#* @apiTitle RUN PROAST v61.5
#* @apiDescription PROAST is a software package that has been developed by the Dutch
#* National Institute for Public Health and the Environment for the statistical ana
#* lysis of dose-response data. It can be used for
#* 1) dose-response modelling,
#* 2) deriving a BenchMark Dose in human risk assessment, and
#* 3) deriving an effect concentration in ecotoxicological risk assessment.
#* For details https://www.rivm.nl/en/proast
# Load training data and parameters
changeV<- read.table('change_variables.txt',col.names='change_variables',colClasses='character')
load('trial1.RData')
library(jsonlite)
library(proast61.5)
# Core predict function
# define parameters with type and description
# name endpoint
# return output as json
#* @post /proast61
#* @json
# original code by GTsiliki (https://github.com/GTsiliki, g_tsiliki@hotmail.com)
enm.proast61 <-
function(dataset,predictionFeature,parameters){
#dataset:= list of 2 objects -
#datasetURI:= character sring, code name of dataset
#dataEntry:= data.frame with 2 columns,
#1st:name of compound,2nd:data.frame with values (colnames are feature names)
#predictionFeature:= character string specifying which is the prediction feature in dataEntry,
#parameters:= list with parameter values-- here awaits a list with four elemnts:
#nTrials (a numeric value indicating number of trials suggested, if 0 then an estimated number is suggested),
#criterion (a character value to indicate which optimal deisgn to apply. Possible values are 'D', 'A', 'I'. Default is 'D'),
#form (a string indicating the formula of the deisgn- 'linear','quad','cubic','cubicS')),
#r2.threshold (numeric value indicating the r2 threshold value. If the data supplied provides r2 value greater
#than the threshold value, a stop message is returned.).
dat<- dataset$dataEntry[,2] # data table
ind.dat<- colnames(dat)
dat1.yname<- predictionFeature # dat1$predictionFeature #string to indicate dependent variable
depend.indx<- which(colnames(dat) %in% dat1.yname)
dat1.xname<- parameters$indiVariable # name
indi.indx<- which(colnames(dat) %in% dat1.xname)
ans.all.new<- ans.all
names.ans.all<- names(ans.all)
pos.ans.all<- which(names(ans.all) %in% changeV[,1])
ans.all.new[[pos.ans.all[1]]]<- indi.indx #xans index
names(ans.all.new[[pos.ans.all[1]]])<- dat1.xname #xans name
ans.all.new[[pos.ans.all[2]]]<- depend.indx #yans index
names(ans.all.new[[pos.ans.all[2]]])<- dat1.yname #yans name
ans.all.new[[pos.ans.all[3]]]<- getwd() #working dir
ans.all.new[[pos.ans.all[4]]]<- dataset$datasetURI #data.name
ans.all.new[[pos.ans.all[5]]]<- dat #odt - dataset
ans.all.new[[pos.ans.all[6]]]<- ncol(dat) #number of columns (i.e. variables) in data set
ans.all.new[[pos.ans.all[7]]]<- colnames(dat) #varnames
ans.all.new[[pos.ans.all[8]]]<- dat1.xname #x name
ans.all.new[[pos.ans.all[9]]]<- dat1.yname #y name
if(range(dat[,depend.indx])[1]<=0){ans.all.new$auto.detlim<- TRUE} #ans.all.new$detlim<- 0.001} # to avoid question for detection limit
ans.all.new$cont <- TRUE #; ans.all<- TRUE
#res0<- f.quick.con(ans.all.new,validate=T)#
#graphics.off()
#.ypos <- 95
devAskNewPage(ask = FALSE)
options(device.ask.default = FALSE)
# run proast & save print out in file
sink("report_from_screen.txt",append=TRUE)
assign(".ypos", 95, envir = .GlobalEnv)
# devAskNewPage(ask = FALSE)
res1<- f.quick.con(ans.all.new,validate=T)#
sink()
#sink.number()
closeAllConnections()
#read print out
rfs<- readLines("report_from_screen.txt")
rfs[which(rfs=="")]<- "\n"
unlink("report_from_screen.txt")
#sink()
#sink.number()
#save, read and resave garphs
g.num<- as.numeric(dev.cur()[1])
c.fig<- as.data.frame(t(character(g.num-1)))
colnames(c.fig)<- c('FittedValues','NestedHillModels','NestedExponentialModels')#paste('figure',1:(g.num-1),sep='')
w<- 0
for(i in g.num:2){
w<-w+1
if(w>3){break}
else{
num1<- g.num-i+1
labeli<- paste('graph',num1,'.png',sep='')
dev.set(which=i)
x<- recordPlot()
png(labeli)
replayPlot(x)
dev.off()
png1<- file(labeli,'rb') #open for reading in binary mode
N<- 1e6
pngfilecontents <- readBin(png1, what="raw", n=N)
close(png1)
#print(c(i))
c.fig[num1]<- fromJSON(toJSON(pngfilecontents))#base64(pngfilecontents)
### NOTE: fromJSON(toJSON( is a round-about since base64 doesn't work with JSON
### and to JSON is doing base64 internally
}
}
if(dev.cur()>1){for(j in 1:(dev.cur()-1)){dev.off()}}
#graphics.off()
liks.mat<- rbind(res1$logliks,res1$EXP$logliks,res1$HILL$logliks)
thres.mat<- c(which(liks.mat[,1]=='--'),dim(liks.mat)[1]+1)
n1.mat<- c('LR',res1$model.sublist)
n2.mat<- character(dim(liks.mat)[1])
w<- 2
while(w<=length(thres.mat)){
n2.mat[thres.mat[w-1]:(thres.mat[w]-1)]<- rep(n1.mat[w-1],thres.mat[w]-thres.mat[w-1])
#n2.len<- thres.mat[w]-thres.mat[w-1]
#n2.mat[thres.mat[w-1]:(thres.mat[w]-1)]<- paste(n2.mat[thres.mat[w-1]:(thres.mat[w]-1)],1:n2.len,sep='')
w<- w+1
}
n3.mat<- paste(n2.mat,liks.mat[,1])
liks.mat.s<- liks.mat
liks.mat.s1<- as.matrix(liks.mat.s)
colnames(liks.mat.s1)<- NULL
x1<- as.list(as.data.frame(t(liks.mat.s1)))
names(x1)<- n3.mat
calc5<- list(colNames=colnames(liks.mat.s),values=x1)
l1<- list(logLik=calc5)
outP<- list(singleCalculations=list(reportFromScreen=unbox(paste(rfs,collapse=' '))),
arrayCalculations=l1,figures=unbox(c.fig))
#unlink('.trial1.RData')
return(outP)
}