使用R從流量數據計算喇叭曲線

Question

我正在尋找一種方法來計算來自靜脈設置流量數據的喇叭曲線（並且最好使用ggplot2，儘管這不重要）。這些是基於來自IV泵的流量數據的曲線。

有關喇叭曲線的示例，請參見http://www.australianprescriber.com/magazine/18/2/49/51/。

任何幫助非常感謝。

皮特

---

下面是一些流數據（速度遠遠超過IV泵 - 從流體變暖泵（I誤所述先前重力流），這是高流動性，但相同類型的圖案）：http://pastebin.com/vJmGcJmn

到目前爲止，我有這樣的：

flow<-read.table(file="flow.dat",header=T) 
flow$diff<-c(0,diff(flow$Mass)) 
ggplot(data=flow, aes(x=Secs,y=diff)) + geom_line() 

注意，所需的小號圖是按照ISO 60601-2-24（如果這意味着什麼！）

該流量設定爲150ml /分鐘。

Answer 1

對不起，這花了這麼久，時間緊迫。

您提供的數據似乎不符合您對小號曲線所表示的內容的描述，或者我錯過了很大的東西。總之，如果你能夠描述需要用數據做什麼，我將不勝感激。

當您設法爲輸出形成數據時，您可以將其粘貼到下面的代碼中，並且它應該會生成一個圖。我會根據您的需求定製它。

# generate some random data 
trump <- data.frame(
    curve1 = rev(sort(rchisq(100, df = 2))) * rnorm(100, mean = 5, sd = 0.1) + 3, 
    curve2 = -rev(sort(rchisq(100, df = 2))) * rnorm(100, mean = 5, sd = 0.1) - 3 
) 
trump <- trump[seq(from = 1, to = 100, by = 3), ] 

ggplot(trump, aes(x = 1:nrow(trump), y = curve1)) + 
    geom_line() + 
    geom_point() + 
    geom_line(aes(y = curve2)) + 
    geom_point(aes(y = curve2)) + 
    geom_hline(aes(yintercept = 0), linetype = "solid") + # overall percentage error 
    geom_hline(aes(yintercept = 5), linetype = "dashed") + # set rate 
    xlab("Observation windows (in minutes)") + 
    ylab("% error") + 
    annotate("text", x = 8, y = -1.5, label = "overall percentage error") + 
    annotate("text", x = 5, y = 3, label = "set rate") + 
    annotate("text", x = 10, y = -24, label = "min. error") + 
    annotate("text", x = 10, y = 24, label = "max. error") + 
    theme_bw() 

Answer 2

非常感謝羅馬對他的幫助。這是我們結束這樣做的方式。

flow<-read.table("iv.txt",comment.char="#",header=TRUE) 
Q<-as.data.frame(c(0,(60*diff(flow$wt,1,1)/0.5))) # Q is used for flow in the standard 
names(Q)<-"Q" 
# Now combine them 
l<-list(flow,Q) 
dflow<-as.data.frame(l) 
t1flow<-subset(dflow, dflow$secs>60*60 & dflow$secs<=120*60) # require the second hour of the data 
t1flow$QE<-((t1flow$Q-setflow)/setflow)*100 # calculate the error 
library(TTR) # use this for moving averages 
# 
# Avert your eyes! I know there must be a slicker way of doing this ..... 
#  
# Calculate the moving average using a sample rate of every 4 readings (2 mins of 30sec readings) 
QE2<-SMA(t1flow$QE,4) 
minQE2<-min(QE2,na.rm=TRUE) 
maxQE2<-max(QE2,na.rm=TRUE) 
# now for the 5 minute window 
QE5<-SMA(t1flow$QE,10) 
minQE5<-min(QE5,na.rm=TRUE) 
maxQE5<-max(QE5,na.rm=TRUE) 
# Set window to 11 mins 
QE11<-SMA(t1flow$QE,22) 
minQE11<-min(QE11,na.rm=TRUE) 
maxQE11<-max(QE11,na.rm=TRUE) 
# Set window to 19 mins 
QE19<-SMA(t1flow$QE,38) 
minQE19<-min(QE19,na.rm=TRUE) 
maxQE19<-max(QE19,na.rm=TRUE) 
# Set window to 31 mins 
QE31<-SMA(t1flow$QE,62) 
minQE31<-min(QE31,na.rm=TRUE) 
maxQE31<-max(QE31,na.rm=TRUE) 
# 
# OK - you can look again :-) 
# 
# create a data frame from this data 
trump<-data.frame(c(2,5,11,19,31),c(minQE2,minQE5,minQE11,minQE19,minQE31),c(maxQE2,maxQE5,maxQE11,maxQE19,maxQE31)) 
names(trump)<-c("T","minE","maxE") 
A<-mean(t1flow$QE) # calculate the overall mean percentage error 
error_caption<-paste("overall percentage error = ",A,"%") # create the string to label the error line 
# plot the graph 
ggplot(trump, aes(x = T, y = minE)) + 
    geom_line() + 
    geom_point(color="red") + 
    geom_line(aes(y = maxE)) + 
    geom_point(aes(y = maxE),colour="red") + 
    geom_hline(aes(yintercept = 0), linetype = "dashed") + # overall percentage error 
    geom_hline(aes(yintercept = A), linetype = "solid") + # set rate 
    xlab("Observation windows (in minutes)") + 
    ylab("% error") + 
    scale_x_continuous(breaks=c(0,2,5,11,19,31),limits=c(0,32)) + # label the x axis only at the window values 
    annotate("text", x = 10, y = A-0.5, label = error_caption) + # add the error line label 
    opts(title="Trumpet curve for Test Data") 

也可hYou可以看到最終圖形[點擊這裏] [2]

Answer 3

我目前正在嘗試解決小號曲線計算錯誤。基於IEC 60601-2-24：1998 - 第八章 - 操作數據的準確性 和防止危險輸出。我在這裏的R代碼裏面：

flow<-read.table("c:\\ciringe.txt",comment.char="#",header=TRUE) 
    #parameters 
    # setflow = 1 
    # P = 1,2,5,11,19,31 

    P1<-1 
    P2<-2 
    P5<-5 
    P11<-11 
    P19<- 19 
    P31<-31 

    P1m = ((60 - P1)/0.5) + 1 
    P2m = ((60 - P2)/0.5) + 1 
    P5m = ((60 - P5)/0.5) + 1 
    P11m = ((60 - P11)/0.5) + 1 
    P19m = ((60 - P19)/0.5) + 1 
    P31m = ((60 - P31)/0.5) + 1 


    setflow<-1 
    mQE1<-0.5 
    mQE2<-0.25 
    mQE5<-0.1 
    mQE11<-0.045 
    mQE19<-0.0263 
    mQE31<-0.0161 

    Q<-as.data.frame(c(0,(60*diff(flow$wt,1,1)/0.5*0.998))) # Q is used for flow in the standard 
    names(Q)<-"Q" 
    # Now combine them 
    l<-list(flow,Q) 
    dflow<-as.data.frame(l) 
    t1flow<-subset(dflow, dflow$secs>=3600 & dflow$secs<=7200) # require the second hour of the data 

    #overall 
    t1flow$QE<-(((t1flow$Q-setflow)/setflow)*100) # calculate the error 

    t1flow$QE1<-(((t1flow$Q-setflow)/setflow)*100) * mQE1 # calculate the error 
    t1flow$QE2<-(((t1flow$Q-setflow)/setflow)*100) * mQE2 # calculate the error 
    t1flow$QE5<-(((t1flow$Q-setflow)/setflow)*100) * mQE5 # calculate the error 
    t1flow$QE11<-(((t1flow$Q-setflow)/setflow)*100) * mQE11 # calculate the error 
    t1flow$QE19<-(((t1flow$Q-setflow)/setflow)*100) * mQE19 # calculate the error 
    t1flow$QE31<-(((t1flow$Q-setflow)/setflow)*100) * mQE31 # calculate the error 

    library(TTR) # use this for moving averages 
    # 
    # Avert your eyes! I know there must be a slicker way of doing this ..... 
    #  
    # Calculate the moving average using a sample rate of every 
    # 4 readings (2 mins of 30sec readings) 

    # now for the 1 minute window 
    QE1<-SMA(t1flow$QE1,2) 
    minQE1<-min(QE1,na.rm=TRUE) 
    maxQE1<-max(QE1,na.rm=TRUE) 

    # now for the 2 minute window 


    QE2<-SMA(t1flow$QE2,4) 
    minQE2<-min(QE2,na.rm=TRUE) 
    maxQE2<-max(QE2,na.rm=TRUE) 

    # now for the 5 minute window 

    QE5<-SMA(t1flow$QE5,10) 
    minQE5<-min(QE5,na.rm=TRUE) 
    maxQE5<-max(QE5,na.rm=TRUE) 

    # Set window to 11 mins 

    QE11<-SMA(t1flow$QE11,22) 
    minQE11<-min(QE11,na.rm=TRUE) 
    maxQE11<-max(QE11,na.rm=TRUE) 

    # Set window to 19 mins 

    QE19<-SMA(t1flow$QE19,38) 
    minQE19<-min(QE19,na.rm=TRUE) 
    maxQE19<-max(QE19,na.rm=TRUE) 

    # Set window to 31 mins 


    QE31<-SMA(t1flow$QE31,62) 
    minQE31<-min(QE31,na.rm=TRUE) 
    maxQE31<-max(QE31,na.rm=TRUE) 

    # 
    # OK - you can look again :-) 
    # 
    # create a data frame from this data 
    trump<- data.frame(c(1,2,5,11,19,31),c(minQE1,minQE2,minQE5,minQE11,minQE19,minQE31), c(maxQE1,maxQE2,maxQE5,maxQE11,maxQE19,maxQE31)) 
    names(trump)<-c("T","minE","maxE") 
    A<-mean(t1flow$QE) # calculate the overall mean percentage error 
    error_caption<-paste("overall percentage error = ",A,"%") # create the string to label the error line 
    # plot the graph 
    library(ggplot2) 
    ggplot(trump, aes(x = T, y = minE)) + 
    geom_line() + 
    geom_point(color="red") + 
    geom_line(aes(y = maxE)) + 
    geom_point(aes(y = maxE),colour="red") + 
    geom_hline(aes(yintercept = 0), linetype = "dashed") + # overall percentage error 
    geom_hline(aes(yintercept = A), linetype = "solid") + # set rate 
    xlab("Observation windows (in minutes)") + 
    ylab("% error") + 
    scale_x_continuous(breaks=c(0,2,5,11,19,31),limits=c(0,32)) + # label the x axis only at the window values 
    annotate("text", x = 10, y = A-0.5, label = error_caption) # add the error line label 

其實，我不明白總和系列在ISO文檔陳述（百分比變化）

Answer 4

你的公式不ISO 60601-2-24匹配方程。該規範有點繁瑣，但沒有2分鐘的窗口或移動平均線。首先，將數據格式化爲1分鐘樣本（根據規範）。做一個1分鐘的間隔和列的數據陣列：

colTotalLiquid = 1 #total liquid (precalculated from weight) 
colProgRate = 2   #prog rate 
colQi = 3   # rate Q(i) from IEC60601-2-24 


obsWindow = (2, 5, 11, 19, 31)     # observation windows in minutes 
analysisT1start = 60       # start of analysis period T1 
analysisT1end = 120        # end of analysis period T1 

t1err = dict.fromkeys(obsWindow, None) 

progRate = data[analysisT1start][colProgRate] 
A = 100 * (((data[analysisT1end][colTotalLiquid] - data[analysisT1start][colTotalLiquid]) * hours/(analysisT1end - analysisT1start)) - progRate)/progRate 

t1TrumpetX = []  #mintes 
t1TrumpetY1 = [] #Ep(max) 
t1TrumpetY2 = [] #Ep(min) 
t1TrumpetY3 = [] #mean err 
t1TrumpetY4 = [] #prog Rate 

for p in obsWindow: 
    m = (analysisT1end - analysisT1start) - p + 1 
    t1err[p] = {'m': m} 
    EpMax = 0 
    EpMin = 0 
    for j in range(1, m + 1): 
     errSum = 0 
     for i in range(j, j + p): 
      errSum += (1.0/p) * 100 * (data[analysisT1start + i][colQi] - progRate)/progRate 
     if errSum > EpMax: 
      EpMax = errSum 
     if errSum < EpMin: 
      EpMin = errSum 
    t1err[p]['EpMax'] = EpMax 
    t1err[p]['EpMin'] = EpMin 
    t1TrumpetX.append(p) 
    t1TrumpetY1.append(EpMax) 
    t1TrumpetY2.append(EpMin) 
    t1TrumpetY3.append(A) 
    t1TrumpetY4.append(0) 

tplot = PdfPages('trumpet curve.pdf') 
p1 = plt.figure() 
plt.plot(t1TrumpetX, t1TrumpetY1) 
plt.plot(t1TrumpetX, t1TrumpetY2) 
plt.plot(t1TrumpetX, t1TrumpetY3) 
plt.plot(t1TrumpetX, t1TrumpetY4, '--') 
plt.legend(('Ep(max)','Ep(min)', 'overall error', 'set rate'), fontsize='xx-small', loc='best') 
plt.title(baseName + ': IEC60601-2-24 Trumpet Curve for Second Hour', fontsize='x-small') 
plt.xlabel('Observation window (min)', fontsize='x-small') 
plt.ylabel('Percentage error of flow', fontsize='x-small') 
plt.ylim((-15, 15)) 
ax = plt.axes() 
ax.set_xticks(obsWindow) 
tplot.savefig(p1) 
plt.close(p1) 
tplot.close(p1) 

這應該得到你的標準的喇叭曲線。沒有必要在整體錯誤A上加上數字標籤，但可以用註釋來完成。