0

I am currently working on a class project where I am supposed to analyze some results in few select papers. I want to calculate the sensitivity indices of the parameters with respect to the basic reproduction number (R0) using the formula S(p)=(p/R0)*(∂R0/∂p) where p is a parameter and S(p) is the sensitivity index of p. I have already calculated the indices manually, but I was wondering if there is a way to automate this process using R. The formula for R0 and the parameter values are given below.

beta_s = 0.274,
alpha_a = 0.4775,
alpha_u = 0.695,
mu = 0.062,
q_i = 0.078,
gamma_a = 0.29,
1/eta_i = 0.009,
1/eta_u = 0.05


R0 = (beta_s*alpha_a)/(gamma_a+mu) + (beta_s*alpha_u*gamma_a*(1-q_i))/((gamma_a+mu)*(eta_u+mu))

It would be great if someone could help me with finding the sensitivity indices using R. Thanks a lot for your time!

Edited code based on @epsilonG's answer

beta_s = 0.274
alpha_a = 0.4775
alpha_u = 0.695
mu = 0.062
q_i = 0.078
gamma_a = 0.29
1/eta_i = 0.009
1/eta_u = 0.05

R0 = (beta_s*alpha_a)/(gamma_a+mu) + (beta_s*alpha_u*gamma_a*(1-q_i))/((gamma_a+mu)*(eta_u+mu))

# Create function
func <- expression((beta_s*alpha_a)/(gamma_a+mu) + (beta_s*alpha_u*gamma_a*(1-q_i))/((gamma_a+mu)*(eta_u+mu))

# Calculate the derivative of R0 to alpha_a

dalpha_a <- deriv(func, 'alpha_a')
val <- eval(dalpha_a)

# Calculate the sensitivity index of alpha_a

S_alpha_a <- (alpha_a/R0)*val
Hew123
  • 63
  • 5

2 Answers2

1
# Create function
func <- expression((beta_s*alpha_a)/(gamma_a+mu) + 
        (beta_s*alpha_u*gamma_a*(1-q_i))/((gamma_a+mu)))

# Calculate the derivative of R0 to beta_s
dbeta_s <- deriv(func , 'beta_s')
val <- eval(dbeta_s)

# Calculate the sensitivity index of beta_s
s_beta_s <- beta_s/R0 * val

You may use for loop do this through all the parameters.

epsilonG
  • 36
  • 3
  • Hey, thanks for the code. It actually worked and I got some values. However, the values given in the paper are (βs, μ, γa, αa, ηu, αu, qi) = (1.000, − 0.6060, −0.18464, 0.22348, −0.34666, 0.77652, −0.06474). None of the values that I got from the code matches these even though the formula and the parameter values are correct. – Hew123 Dec 04 '22 at 21:05
0

I figured the issue with this code. In fact, it was just a minor thing, I used the function D, instead of deriv.

beta_s = 0.274
alpha_a = 0.4775
alpha_u = 0.695
mu = 0.062
q_i = 0.078
gamma_a = 0.29
1/eta_i = 0.009
1/eta_u = 0.05

R0 = (beta_s*alpha_a)/(gamma_a+mu) + (beta_s*alpha_u*gamma_a*(1-q_i))/((gamma_a+mu)*(eta_u+mu))

# Create function
func <- expression((beta_s*alpha_a)/(gamma_a+mu) + (beta_s*alpha_u*gamma_a*(1-q_i))/((gamma_a+mu)*(eta_u+mu))

# Calculate the derivative of R0 to alpha_a

dalpha_a <- D(func, 'alpha_a')
val <- eval(dalpha_a)

# Calculate the sensitivity index of alpha_a

S_alpha_a <- (alpha_a/R0)*val
Hew123
  • 63
  • 5