This repository has been archived by the owner on Mar 30, 2020. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ComputeAndSaveWNFit.R
107 lines (104 loc) · 3.52 KB
/
ComputeAndSaveWNFit.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
#' @title Compute and save Wagner-Nelson fit
#'
#' @description Reads record with given BreathTestRecordID from database and computes
#' Wagner-Nelsoen fit parameter which are written back to table BreathTestParameter
#' Existing parameters are overwritten.
#' See Sanaka, Yamamoto, Ishii, Kuyama, Digestion 2004 (69), 71-78
#'
#' @param con connection to SQlite database
#' @param BreathTestRecordID these data will be read from BreathTestTimeSeries
#' @examples
#' sqliteFile = CreateSimulatedBreathTestDatabase()
#' con = OpenSqliteConnection(sqliteFile)
#' BreathTestRecordID = 1
#' ComputeAndSaveWNFit(con,BreathTestRecordID)
#' dbDisconnect(con)
#'
#' @author Dieter Menne, \email{dieter.menne@@menne-biomed.de}
#' @export
#'
ComputeAndSaveWNFit = function(con,BreathTestRecordID) {
k = dbGetQuery(
con, paste0(
"SELECT Value as k from BreathTestParameter where BreathTestRecordID=",
BreathTestRecordID," and Parameter ='k'"
)
)[1,1]
if (is.na(k))
k = 0.65 / 60 # Default value is 0.65, but we use minutes everywhere
data = dbGetQuery(
con,paste0(
"SELECT Time, Value as PDR from BreathTestTimeSeries where BreathTestRecordID=" ,
BreathTestRecordID, " and Parameter = 'PDR' and Time > 0 order by Time"
)
)
if (inherits("data","try-error"))
stop(paste("No PDR data found for BreathTestRecordID",BreathTestRecordID))
auct = AUCt(data$Time,data$PDR)
aucInf = auct[length(auct)] + data$PDR[nrow(data)] / k # equation 7 in Sanaka
frac = 1 - (auct + data$PDR / k) / aucInf
# compute 0.5-crossing
spl = smooth.spline(data$Time,frac)
t50 = try(uniroot(function(t)
predict(spl,t)$y - 0.5,c(0,max(data$Time)))$root,silent = TRUE)
if (class(t50) == "try-error")
t50 = 0
ret = dbExecute(
con, paste0(
"DELETE FROM BreathTestParameter where BreathTestRecordID=",
BreathTestRecordID ," and Method ='WN'"
)
)
# Write parameters and coefficients
pars = data.frame(
BreathTestRecordID = BreathTestRecordID,
Parameter = "t50",
Method = "WN",
Value = t50
)
# *** Todo: remove dirty trick for invalid WN
pars = pars[pars$Value != 0,]
if (nrow(pars) == 0)
return(0) # Do not report error, we can live without
ret = try(
SaveBreathTestParameters(con, pars),
silent = TRUE
)
if (inherits(ret,"try-error"))
stop(paste0("Could not write WN fit parameters for Record ",BreathTestRecordID))
# Compute predicted WN fit
Time = seq(0,max(data$Time), by = 5)
wn = data.frame(
BreathTestTimeSeriesID = as.integer(NA),
BreathTestRecordID = as.character(BreathTestRecordID),
Time = Time,
Parameter = as.character("WN"),
Value = predict(spl,Time)$y,
stringsAsFactors = FALSE
)
# Delete old values
ret = dbExecute(
con, paste0(
"DELETE FROM BreathTestTimeSeries where BreathTestRecordID=",
BreathTestRecordID ," and Parameter ='WN'"
)
)
# Retrieve column names to get the order right
flds = dbListFields(con,"BreathTestTimeSeries")
q = paste0("INSERT INTO BreathTestTimeSeries VALUES($",
paste0(flds, collapse = ",$"),")")
ret = try(
dbExecute(con, q, params = wn), silent = TRUE
)
if (inherits(ret, "try-error"))
stop(paste0(
"Could not write predicted Wagner-Nelson data for record ",
BreathTestRecordID, "\n", as.character(ret)
))
t50
}
AUCt = function(x, y)
{
# http://stackoverflow.com/questions/14201375/how-can-i-calculate-the-95-credible-limit-of-an-area-under-a-curve-in-r/14204138#14204138
cumsum(c(0,diff(x) * (y[-1] + y[-length(y)])) / 2)
}