Calculating sensitivity indices using R

Question

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

Answer 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.

Answer 2

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