我尝试过使用optim（）优化模型参数，但得到的参数相同。如何获得优化的参数？

library(deSolve) library(reshape2) library(ggplot2) library(ggpubr) # Model Input N <- 1380004385 p <- .87 #0.88689 initial_state_values <- c(S = N - 1, E = 0, A = 0, I = 1, T = 0, F = 0, R = 0) parameters <- c(pie = .6782, mu = 0.0000391, beta =.89768, alpha = 0.24757, phi = 0.08, h = 1.08, mu_I = 0.06891, mu_T = 0.06891, gamma_I = 0.05090, gamma_T = 0.07048) #67446.82054) time <- seq(1, 164, .1) # Input function: Differential equation SEAITFR_fn <- function(time, initial_state_values, parameters) { with(as.list(c(initial_state_values, parameters)), { N <- S+E+A+I+T+F+R lamda <- beta/N * (A + I) dS <- pie - (lamda + mu) * S dE <- lamda * S - (alpha + mu) * E dA <- alpha * (1-p) * E - (mu + phi) * A dI <- alpha * p * E + phi * A -(h + gamma_I + mu_I + mu) * I dT <- h * I - (gamma_T + mu_T) * T dF <- mu_I * I + mu_T * T dR <- gamma_I * I + gamma_T * T - mu * R return(list(c(dS, dE, dA, dI, dT, dF, dR))) }) } output <- as.data.frame(ode(y = initial_state_values, times = time, func = SEAITFR_fn, parms = parameters) ) #output output$total_prevalence <- output$I + output$T #output$total_prevalence # Distance Function SEAITFR_SSQ <- function(parameters, data) { output <- as.data.frame(ode(y = initial_state_values, # vector of initial state times = time, # vector of times func = SEAITFR_fn, parms = parameters) ) data <- na.omit(data) deltas_square <- (output$total_prevalence[output$time %in% data$time] - data$new_cases)^2 SSQ <- sum(deltas_square) return(SSQ) } # Real world data covid_19_data <- read_excel("covid.xlsx") covid_19_data #Optimization optimised <- optim(par = c(pie =.6782, mu = 0.0000391, beta = 0.88689, alpha = 0.24757, phi = 0.08, h = 1.08, mu_I = 0.06891, mu_T = 0.06891, gamma_I = 0.05090, gamma_T = 0.07048), # these are the starting beta # and gamma that will be fed # first, into SSQ_fn fn = SEAITFR_SSQ, data = covid_19_data # this argument comes under "..." # "Further arguments to be passed to fn a ) optimised optimised_model <- as.data.frame(ode(y = initial_state_values, times = time, func = SEAITFR_fn, parms = optimised$par)) # Optimised_model optimised_model$prevalence <- optimised_model$I + optimised_model$T #optimised_model$prevalence # Plotting plot2 <- ggplot() + geom_line(data = optimised_model, aes(x = time, y = prevalence)) + geom_point(data = covid_19_data, aes(x = time, y = new_cases), colour = "red" ) + xlab("Times(days)") + ylab("No. of infection") + labs(title = paste("Calibration of SIR model with optimised value ")) + ylim(0,30000) plot2

2条回答

网友

1楼 · 编辑于 2024-04-25 05:48:16

我能够根据输出让您的模型运行different parameters。此外，我还使用了maximum likelihood estimation，并且它与covid data配合得很好

网友

2楼 · 编辑于 2024-04-25 05:48:16

require(deSolve)
SIR_fn <- function(time, state, parameters) { +

```
with(as.list(c(state, parameters)), {
```
```
    N  <- S+I+R
```
```
    dS <- -beta*S*I/N
```
```
    dI <- beta*S*I/N-gamma*I
```
```
    dR <- gamma*I
```
```
    return(list(c(dS, dI, dR)))
```
```
})
```
}

SIR_SSQ <- function(parameters, dat) { # parameters must contain beta & gamma

# calculate model output using your SIR function with ode()

result <- as.data.frame(ode(y = initial_state_values  # vector of initial state

                                                      # values, with named elements

                          , times = times             # vector of times

                          , func = SIR_fn             # your predefined SIR function

                          , parms = parameters)       # the parameters argument

                                                      # entered with SIR_SSQ()

```
)
```

# SSQ calculation: needs the dat argument (the observed data you are fitting to)

# assumes the data you are fitting to has a column "I"

# select only complete cases, i.e. rows with no NAs, from the dataframe

```
dat <- na.omit(dat)  
```

# select elements where results$time is in dat$time

deltas2 <- (result$I[result$time %in% dat$time] - dat$I)^2

```
SSQ   <- sum(deltas2)
```
```
return(SSQ)
```
}

flu_dat <- read.csv("flu.csv")
head(flu_dat,10) time I 1 1 16928 2 2 21389 3 3 20006 4 4 21368 5 5 25349 6 6 21405 7 7 17662 8 8 17492 9 9 18361 10 10 20951
initial_state_values <- c(S = 100000, I = 1, R = 0)
choose values to start your optimisation
beta_start <- 0.1 gamma_start <- 0.1
times - dense timesteps for a more detailed solution
times <- seq(from = 0, to = 107, by = 0.01)
optim
you will need to run the cells to assign your functions first
optimised <- optim(par = c(beta = beta_start

                    , gamma = gamma_start)      # these are the starting beta

                                                # and gamma that will be fed

                                                # first, into SSQ_fn

```
                    , fn = SIR_SSQ
```

                    , dat = flu_dat  # this argument comes under "..."

                                     # "Further arguments to be passed to fn and gr"

```
  )
```

optimised #have a look at the model output $par beta gamma 0.40487383 0.01637813

美元价值 [1] 25002879441

美元计数函数梯度 81纳

$convergence [1] 0

$message 空的

optimised$par beta gamma 0.40487383 0.01637813
opt_mod <- as.data.frame(ode(y = initial_state_values # named vector of initial

                                                   # state values

                        , times = times            # vector of times

                        ,  func = SIR_fn           # your predefined SIR function

                        , parms = optimised$par))

plot your optimised model output, with the epidemic data using ggplot
require(ggplot2)
opt_plot <- ggplot() opt_plot <- opt_plot + geom_point(aes(x = time, y = I)

                            , colour = "red"

```
                            ,size = 3
```

                            , shape  = "x"

                            , data = flu_dat)

opt_plot <- opt_plot + geom_line(aes(x = time, y = I)

                             , colour = "blue"

                             , data   = opt_mod)

opt_plot
flu_dat time I 1 1 16928 2 2 21389 3 3 20006 4 4 21368 5 5 25349 6 6 21405 7 7 17662 8 8 17492 9 9 18361 10 10 20951 11 11 23042 12 12 24914 13 13 25268 14 14 19019 15 15 19456 16 16 23693 17 17 26499 18 18 27867 19 19 30981 20 20 31917 21 21 25968 22 22 30676 23 23 36208 24 24 34426 25 25 40356 26 26 45353 27 27 41359 28 28 40409 29 29 40057 30 30 45980 31 31 51488 32 32 55588 33 33 51841 34 34 45549 35 35 49580 36 36 44162 37 37 60755 38 38 59740 39 39 62567 40 40 67859 41 41 60038 42 42 59085 43 43 58793 44 44 67427 45 45 67635 46 46 77102 47 47 71672 48 48 62470 49 49 60716 50 50 61477 51 51 64239 52 52 71915 53 53 68511 54 54 73208 55 55 65375 56 56 54864 57 57 55947 58 58 65785 59 59 71832 60 60 67599 61 61 67911 62 62 57839 63 63 46269 64 64 44567 65 65 57283 66 66 54433 67 67 59481 68 68 58322 69 69 54449 70 70 46369 71 71 48733 72 72 46910 73 73 56811 74 74 51742 75 75 64764 76 76 46790 77 77 41181 78 78 36484 79 79 45034 80 80 47393 81 81 44087 82 82 48253 83 83 43653 84 84 34510 85 85 36432 86 86 39979 87 87 45520 88 88 45120 89 89 46929 90 90 45414 91 91 34787 92 92 35138 93 93 41677 94 94 40951 95 95 43781 96 96 50129 97 97 43192 98 98 31511 99 99 23545 100 100 26845 101 101 33895 102 102 36066 103 103 47439 104 104 41003 105 105 34305 106 106 33842 107 107 39385

choose values to start your optimisation

times - dense timesteps for a more detailed solution

optim

you will need to run the cells to assign your functions first

plot your optimised model output, with the epidemic data using ggplot

相关问题更多 >

编程相关推荐

热门问题

热门文章