0

我使用包 deSolve 在 R 中制作了一个 ODE 模型。目前，该模型的输出为我提供了“观察到的”疾病流行率（即流行率不考虑诊断缺陷）。

但是，我想调整模型以输出“真实”流行率，使用称为 Rogan-Gladen 估计器（http://influentialpoints.com/Training/estimating_true_prevalence.htm）的简单调整公式：

真实患病率 = (表观上一个 + (Specificity-1)) / (Specificity + (Sensitivity-1))

正如您将在下面的代码中看到的那样，我只尝试调整其中一个微分方程 (diggP)。

在没有调整的情况下运行模型会给出预期的输出（0 到 1 之间的比例）。但是，尝试使用 RG 估计器调整模型会产生虚假输出（比例小于 0）。

任何关于这里可能出现问题的建议将不胜感激。

# Load required packages
library(tidyverse)
library(broom)
library(deSolve)

# Set time (age) for function
time = 1:80

# Defining exponential decay of lambda over age
y1 = 0.003 + (0.15 - 0.003) * exp(-0.05 * time) %>% jitter(10)
df <- data.frame(t = time, y = y1)
fit <- nls(y ~ SSasymp(time, yf, y0, log_alpha), data = df)
fit

# Values of lambda over ages 1-80 years
data <- as.matrix(0.003 + (0.15 - 0.003) * exp(-0.05 * time))

lambda<-as.vector(data[,1])
t<-as.vector(seq(1, 80, by=1))
foi<-cbind(t, lambda)
foi[,1]

# Making lambda varying by time useable in the ODE model
input <- approxfun(x = foi[,1], y = foi[,2], method = "constant", rule = 2)


# Model
ab <- function(time, state, parms) {
  with(as.list(c(state, parms)), {

    # lambda, changing by time
    import<-input(time) 

    # Derivatives

    # RG estimator: 
    #True prevalence = (apparent prev + (sp-1)) / (sp + (se-1))

    diggP<- (((import * iggN) - iggR * iggP) + (sp_igg-1)) / (sp_igg + (se_igg-1))
    diggN<- (-import*iggN) + iggR*iggP

    dtgerpP<- (0.5*import)*tgerpN -tgerpR*tgerpP
    dtgerpN<- (0.5*-import)*tgerpN + tgerpR*tgerpP

    # Return results
    return(list(c(diggP, diggN, dtgerpP, dtgerpN))) 
  })
}

# Initial values
yini  <- c(iggP=0, iggN=1, 
           tgerpP=0, tgerpN=1) 

# Parameters
pars  <- c(iggR = 0, tgerpR = (1/8)/12, 
           se_igg = 0.95, sp_igg = 0.92)

# Solve model
results<- ode(y=yini, times=time, func=ab, parms = pars)

# Plot results
plot(results, xlab="Time (years)", ylab="Proportion")
4

0 回答 0