Настройка вывода модели ODE с использованием оценки Рогана-Глейдена в R - PullRequest
Новая
       71

Настройка вывода модели ODE с использованием оценки Рогана-Глейдена в R

0 голосов
/ 01 апреля 2019

Я сделал модель ODE в R, используя пакет deSolve.В настоящее время выходные данные модели дают мне «наблюдаемую» распространенность заболевания (т.е. распространенность, не учитывающую диагностическое несовершенство).

Однако я хочу настроить модель для вывода «истинной» распространенности, используяпростая формула корректировки, называемая оценщиком Рогана-Глейдена (http://influentialpoints.com/Training/estimating_true_prevalence.htm):

Истинная распространенность = (Кажущаяся прев. + (Специфичность-1)) / (Специфичность + (Чувствительность-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")
...