Hello All,
I am currently looking on a transmission model for STD transmission within a
population. I am able to run my full code and the ODE function, but when I
look at my output, all I get is "NA" for each time step beyond the first.
There doesn't seem to be any syntax error, and I do get my entire program to
run. Here is my code:
setwd("C:/Users/L/Documents/MastersThesis")
require(deSolve);
########
#Model 1
########
#The function
HPVInfection<-function(t,y,p){
XFL = y[1]; #number of susceptible unvaccinated females low risk
XFM = y[2]; #number of susceptible unvaccinated females medium risk
XFH = y[3]; #number of susceptible unvaccinated females high risk
XML = y[4]; #number of susceptible unvaccinated males low risk
XMM = y[5]; #number of susceptible unvaccinated males medium risk
XMH = y[6]; #number of susceptible unvaccinated males high risk
Y1FL = y[7]; #number of infected unvaccinated females low risk infected
with vaccine strain
Y1FM = y[8]; #number of infected unvaccinated females medium risk low
risk
infected with vaccine strain
Y1FH = y[9]; #number of infected unvaccinated females high risk low risk
infected with vaccine strain
Y1ML = y[10]; #number of infected unvaccinated males low risk low risk
infected with vaccine strain
Y1MM = y[11]; #number of infected unvaccinated males medium risk low
risk
infected with vaccine strain
Y1MH = y[12]; #number of infected unvaccinated males high risk low risk
infected with vaccine strain
Y2FL = y[13]; #number of infected unvaccinated females low risk infected
with non-vaccine strain
Y2FM = y[14]; #number of infected unvaccinated females medium risk low
risk
infected with non-vaccine strain
Y2FH = y[15]; #number of infected unvaccinated females high risk low
risk
infected with non-vaccine strain
Y2ML = y[16]; #number of infected unvaccinated males low risk low risk
infected with non-vaccine strain
Y2MM = y[17]; #number of infected unvaccinated males medium risk low
risk
infected with non-vaccine strain
Y2MH = y[18]; #number of infected unvaccinated males high risk low risk
infected with non-vaccine strain
ZFL = y[19]; #number of immune females low risk
ZFM = y[20]; #number of immune females medium risk
ZFH = y[21]; #number of immune females high risk
ZML = y[22]; #number of immune males low risk
ZMM = y[23]; #number of immune males medium risk
ZMH = y[24]; #number of immune males high risk
VFL = y[25]; #number of susceptible vaccinated females low risk
VFM = y[26]; #number of susceptible vaccinated females medium risk
VFH = y[27]; #number of susceptible vaccinated females high risk
VML = y[28]; #number of susceptible vaccinated males low risk
VMM = y[29]; #number of susceptible vaccinated males medium risk
VMH = y[30]; #number of susceptible vaccinated males high risk
W1FL = y[31]; #number of infected vaccinated females low risk infected
with
vaccine strain
W1FM = y[32]; #number of infected vaccinated females medium risk
infected
with vaccine strain
W1FH = y[33]; #number of infected vaccinated females high risk infected
with vaccine strain
W1ML = y[34]; #number of infected vaccinated males low risk infected
with
vaccine strain
W1MM = y[35]; #number of infected vaccinated males medium risk infected
with vaccine strain
W1MH = y[36]; #number of infected vaccinated males high risk infected
with
vaccine strain
W2FL = y[37]; #number of infected vaccinated females low risk infected
with
non-vaccine strain
W2FM = y[39]; #number of infected vaccinated females medium risk
infected
with non-vaccine strain
W2FH = y[40]; #number of infected vaccinated females high risk infected
with non-vaccine strain
W2ML = y[41]; #number of infected vaccinated males low risk infected
with
non-vaccine strain
W2MM = y[42]; #number of infected vaccinated males medium risk infected
with non-vaccine strain
W2MH = y[43]; #number of infected vaccinated males high risk infected
with
non-vaccine strain
with(as.list(p), {
dXFL.dt = (0.5 * mew * omega[1,1] * (1-phi) * total) -
((partner[1,1] *
beta[1,1] * ((((Y1ML + Y2ML + (tau[1,1] * W1ML) + (tau[1,2] * W2ML)) /
population[1,2]) * rho[1,1]) + (((Y1MM + Y2MM + (tau[1,1]*W1MM) + (tau[1,2]
* W2MM))/population[2,2]) * rho[1,2]) + (((Y1MH + Y2MH + (tau[1,2] * W1MH) +
(tau[1,2] * W2MH)) / population[3,2]) * rho[1,3])) + mew) * XFL) + (sigma *
VFL);
dXFM.dt = (0.5 * mew * omega[2,1] * (1-phi) * total) -
((partner[2,1] *
beta[1,1] * ((((Y1ML + Y2ML + (tau[1,1] * W1ML) + (tau[1,2] * W2ML)) /
population[1,2]) * rho[2,1]) + (((Y1MM + Y2MM + (tau[1,1]*W1MM) + (tau[1,2]
* W2MM))/population[2,2]) * rho[2,2]) + (((Y1MH + Y2MH + (tau[1,1] * W1MH) +
(tau[1,2] * W2MH)) / population[3,2]) * rho[2,3])) + mew) * XFM) + (sigma *
VFM);
dXFH.dt = (0.5 * mew * omega[3,1] * (1-phi) * total) -
((partner[3,1] *
beta[1,1] * ((((Y1ML + Y2ML + (tau[1,1] * W1ML) + (tau[1,2] * W2ML)) /
population[1,2]) * rho[3,1]) + (((Y1MM + Y2MM + (tau[1,1]*W1MM) + (tau[1,2]
* W2MM))/population[2,2]) * rho[3,2]) + (((Y1MH + Y2MH + (tau[1,1] * W1MH) +
(tau[1,2] * W2MH)) / population[3,2]) * rho[3,3])) + mew) * XFH) + (sigma *
VFH);
dXML.dt = (0.5 * mew * omega[1,1] * (1-phi) * total) -
((partner[1,1] *
beta[2,1] * ((((Y1FL + Y2FL + (tau[2,1] * W1FL) + (tau[2,2] * W2FL)) /
population[1,1]) * rho[1,1]) + (((Y1FM + Y2FM + (tau[2,1]*W1FM) + (tau[2,2]
* W2FM))/population[2,1]) * rho[1,2]) + (((Y1FH + Y2FH + (tau[2,1] * W1FH) +
(tau[2,2] * W2FH)) / population[3,1]) * rho[1,3])) + mew) * XML) + (sigma *
VML);
dXMM.dt = (0.5 * mew * omega[2,1] * (1-phi) * total) -
((partner[2,1] *
beta[2,1] * ((((Y1FL + Y2FL + (tau[2,1] * W1FL) + (tau[2,2] * W2FL)) /
population[1,1]) * rho[2,1]) + (((Y1FM + Y2FM + (tau[2,1]*W1FM) + (tau[2,2]
* W2FM))/population[2,1]) * rho[2,2]) + (((Y1FH + Y2FH + (tau[2,1] * W1FH) +
(tau[2,2] * W2FH)) / population[3,1]) * rho[2,3])) + mew) * XMM) + (sigma *
VMM);
dXMH.dt = (0.5 * mew * omega[3,1] * (1-phi) * total) -
((partner[3,1] *
beta[2,1] * ((((Y1FL + Y2FL + (tau[2,1] * W1FL) + (tau[2,2] * W2FL)) /
population[1,1]) * rho[3,1]) + (((Y1FM + Y2FM + (tau[2,1]*W1FM) + (tau[2,2]
* W2FM))/population[2,1]) * rho[3,2]) + (((Y1FH + Y2FH + (tau[2,1] * W1FH) +
(tau[2,2] * W2FH)) / population[3,1]) * rho[3,3])) + mew) * XMH) + (sigma *
VMH);
dY1FL.dt = (XFL * (partner[1,1] * beta[1,1] * ((rho[1,1] *
((Y1ML +
(tau[1,1] * W1ML)) / population[1,2])) + (rho[1,2] * ((Y1MM + (tau[1,1] *
W1MM)) / population[2,2])) + (rho[1,3] * ((Y1MH + (tau[1,1] * W1MH)) /
population[3,2]))))) - ((mew + gamma[2,1]) * Y1FL);
dY1FM.dt = (XFM * (partner[2,1] * beta[1,1] * ((rho[2,1] *
((Y1ML +
(tau[1,1] * W1ML)) / population[1,2])) + (rho[2,2] * ((Y1MM + (tau[1,1] *
W1MM)) / population[2,2])) + (rho[2,3] * ((Y1MH + (tau[1,1] * W1MH)) /
population[3,2]))))) - ((mew + gamma[2,1]) * Y1FM);
dY1FH.dt = (XFH * (partner[3,1] * beta[1,1] * ((rho[3,1] *
((Y1ML +
(tau[1,1] * W1ML)) / population[1,2])) + (rho[2,2] * ((Y1MM + (tau[1,1] *
W1MM)) / population[2,2])) + (rho[3,3] * ((Y1MH + (tau[1,1] * W1MH)) /
population[3,2]))))) - ((mew + gamma[2,1]) * Y1FH);
dY1ML.dt = (XML * (partner[1,1] * beta[2,1] * ((rho[1,1] *
((Y1FL +
(tau[2,1] * W1FL)) / population[1,1])) + (rho[1,2] * ((Y1FM + (tau[2,1] *
W1FM)) / population[2,1])) + (rho[1,3] * ((Y1FH + (tau[2,1] * W1FH)) /
population[3,1]))))) - ((mew + gamma[1,1]) * Y1MM);
dY1MM.dt = (XMM * (partner[2,1] * beta[2,1] * ((rho[2,1] *
((Y1FL +
(tau[2,1] * W1FL)) / population[1,1])) + (rho[2,2] * ((Y1FM + (tau[2,1] *
W1FM)) / population[2,1])) + (rho[2,3] * ((Y1FH + (tau[2,1] * W1FH)) /
population[3,1]))))) - ((mew + gamma[1,1]) * Y1MM);
dY1MH.dt = (XMH * (partner[3,1] * beta[2,1] * ((rho[3,1] *
((Y1FL +
(tau[2,1] * W1FL)) / population[1,1])) + (rho[2,2] * ((Y1FM + (tau[2,1] *
W1FM)) / population[2,1])) + (rho[3,3] * ((Y1FH + (tau[2,1] * W1MH)) /
population[3,1]))))) - ((mew + gamma[1,1]) * Y1MH);
dY2FL.dt = (XFL * (partner[1,1] * beta[1,1] * ((rho[1,1] *
((Y2ML +
(tau[1,2] * W2ML)) / population[1,2])) + (rho[1,2] * ((Y1MM + (tau[1,2] *
W2MM)) / population[2,2])) + (rho[1,3] * ((Y2MH + (tau[1,2] * W2MH)) /
population[3,2]))))) - ((mew + gamma[2,2]) * Y1FL);
dY2FM.dt = (XFM * (partner[2,1] * beta[1,1] * ((rho[2,1] *
((Y2ML +
(tau[1,2] * W2ML)) / population[1,2])) + (rho[2,2] * ((Y1MM + (tau[1,2] *
W2MM)) / population[2,2])) + (rho[2,3] * ((Y2MH + (tau[1,2] * W2MH)) /
population[3,2]))))) - ((mew + gamma[2,2]) * Y1FM);
dY2FH.dt = (XFH * (partner[3,1] * beta[1,1] * ((rho[3,1] *
((Y2ML +
(tau[1,2] * W2ML)) / population[1,2])) + (rho[2,2] * ((Y1MM + (tau[1,2] *
W2MM)) / population[2,2])) + (rho[3,3] * ((Y2MH + (tau[1,2] * W2MH)) /
population[3,2]))))) - ((mew + gamma[2,2]) * Y1FH);
dY2ML.dt = (XML * (partner[1,1] * beta[2,1] * ((rho[1,1] *
((Y2FL +
(tau[2,2] * W2FL)) / population[1,1])) + (rho[1,2] * ((Y1FM + (tau[2,2] *
W2FM)) / population[2,1])) + (rho[1,3] * ((Y2FH + (tau[2,2] * W2FH)) /
population[3,1]))))) - ((mew + gamma[1,2]) * Y1MM);
dY2MM.dt = (XMM * (partner[2,1] * beta[2,1] * ((rho[2,1] *
((Y2FL +
(tau[2,2] * W2FL)) / population[1,1])) + (rho[2,2] * ((Y1FM + (tau[2,2] *
W2FM)) / population[2,1])) + (rho[2,3] * ((Y2FH + (tau[2,2] * W2FH)) /
population[3,1]))))) - ((mew + gamma[1,2]) * Y1MM);
dY2MH.dt = (XMH * (partner[3,1] * beta[2,1] * ((rho[3,1] *
((Y2FL +
(tau[2,2] * W2FL)) / population[1,1])) + (rho[2,2] * ((Y1FM + (tau[2,2] *
W2FM)) / population[2,1])) + (rho[3,3] * ((Y2FH + (tau[2,2] * W2MH)) /
population[3,1]))))) - ((mew + gamma[1,2]) * Y1MH);
dZFL.dt = ((gamma[2,1] * Y1FL) + (gamma[2,2] * Y2FL) +
(gamma[4,1] * W1FL)
+ (gamma[4,2]*W2FL)) - (mew * ZFL);
dZFM.dt = ((gamma[2,1] * Y1FM) + (gamma[2,2] * Y2FM) +
(gamma[4,1] * W1FM)
+ (gamma[4,2]*W2FM)) - (mew * ZFM);
dZFH.dt = ((gamma[2,1] * Y1FH) + (gamma[2,2] * Y2FH) +
(gamma[4,1] * W1FH)
+ (gamma[4,2]*W2FH)) - (mew * ZFH);
dZML.dt = ((gamma[1,1] * Y1ML) + (gamma[1,2] * Y2ML) +
(gamma[3,1] * W1FM)
+ (gamma[3,2]*W2ML)) - (mew * ZML);
dZMM.dt = ((gamma[1,1] * Y1MM) + (gamma[1,2] * Y2MM) +
(gamma[3,1] * W1MM)
+ (gamma[3,2]*W2MM)) - (mew * ZMM);
dZMH.dt = ((gamma[1,1] * Y1MH) + (gamma[1,2] * Y2MH) +
(gamma[3,1] * W1MH)
+ (gamma[3,2]*W2MH)) - (mew * ZMH);
dVFL.dt = (0.5 * mew * omega[1,1] * phi * total) -
((partner[1,1] *
beta[1,1] * ((((Y1ML + Y2ML + (delta[1,1] * tau[1,1] * W1ML) + (delta[1,2] *
tau[1,2] * W2ML)) / population[1,2]) * rho[1,1]) + (((Y1MM + Y2MM +
(delta[1,1] * tau[1,1]*W1MM) + (delta[1,2] * tau[1,2] *
W2MM))/population[2,2]) * rho[1,2]) + (((Y1MH + Y2MH + (delta[1,1] *
tau[1,2] * W1MH) + (delta[1,2] * tau[1,2] * W2MH)) / population[3,2]) *
rho[1,3])) + mew) * VFL) - (sigma * VFL);
dVFM.dt = (0.5 * mew * omega[2,1] * phi * total) -
((partner[2,1] *
beta[1,1] * ((((Y1ML + Y2ML + (delta[1,1] * tau[1,1] * W1ML) + (delta[1,2] *
tau[1,2] * W2ML)) / population[1,2]) * rho[2,1]) + (((Y1MM + Y2MM +
(delta[1,1] * tau[1,1]*W1MM) + (delta[1,2] * tau[1,2] *
W2MM))/population[2,2]) * rho[2,2]) + (((Y1MH + Y2MH + (delta[1,1] *
tau[1,1] * W1MH) + (delta[1,2] * tau[1,2] * W2MH)) / population[3,2]) *
rho[2,3])) + mew) * VFM) - (sigma * VFM);
dVFH.dt = (0.5 * mew * omega[3,1] * phi * total) -
((partner[3,1] *
beta[1,1] * ((((Y1ML + Y2ML + (delta[1,1] * tau[1,1] * W1ML) + (delta[1,2] *
tau[1,2] * W2ML)) / population[1,2]) * rho[3,1]) + (((Y1MM + Y2MM +
(delta[1,1] * tau[1,1]*W1MM) + (delta[1,2] * tau[1,2] *
W2MM))/population[2,2]) * rho[3,2]) + (((Y1MH + Y2MH + (delta[1,1] *
tau[1,1] * W1MH) + (delta[1,2] * tau[1,2] * W2MH)) / population[3,2]) *
rho[3,3])) + mew) * VFH) - (sigma * VFH);
dVML.dt = (0.5 * mew * omega[1,1] * phi * total) -
((partner[1,1] *
beta[2,1] * ((((Y1FL + Y2FL + (delta[2,1] * tau[2,1] * W1FL) + (delta[2,2] *
tau[2,2] * W2FL)) / population[1,1]) * rho[1,1]) + (((Y1FM + Y2FM +
(delta[2,1] * tau[2,1]*W1FM) + (delta[2,2] * tau[2,2] *
W2FM))/population[2,1]) * rho[1,2]) + (((Y1FH + Y2FH + (delta[2,1] *
tau[2,1] * W1FH) + (delta[2,2] * tau[2,2] * W2FH)) / population[3,1]) *
rho[1,3])) + mew) * VML) - (sigma * VML);
dVMM.dt = (0.5 * mew * omega[2,1] * phi * total) -
((partner[2,1] *
beta[2,1] * ((((Y1FL + Y2FL + (delta[2,1] * tau[2,1] * W1FL) + (delta[2,2] *
tau[2,2] * W2FL)) / population[1,1]) * rho[2,1]) + (((Y1FM + Y2FM +
(delta[2,1] * tau[2,1]*W1FM) + (delta[2,2] * tau[2,2] *
W2FM))/population[2,1]) * rho[2,2]) + (((Y1FH + Y2FH + (delta[2,1] *
tau[2,1] * W1FH) + (delta[2,2] * tau[2,2] * W2FH)) / population[3,1]) *
rho[2,3])) + mew) * VMM) - (sigma * VMM);
dVMH.dt = (0.5 * mew * omega[3,1] * phi * total) -
((partner[3,1] *
beta[2,1] * ((((Y1FL + Y2FL + (delta[2,1] * tau[2,1] * W1FL) + (delta[2,2] *
tau[2,2] * W2FL)) / population[1,1]) * rho[3,1]) + (((Y1FM + Y2FM +
(delta[2,1] * tau[2,1]*W1FM) + (delta[2,2] * tau[2,2] *
W2FM))/population[2,1]) * rho[3,2]) + (((Y1FH + Y2FH + (delta[2,1] *
tau[2,1] * W1FH) + (delta[2,2] * tau[2,2] * W2FH)) / population[3,1]) *
rho[3,3])) + mew) * VMH) - (sigma * VMH);
dW1FL.dt = (VFL * (partner[1,1] * beta[1,1] * ((rho[1,1] *
((Y1ML +
(delta[1,1] * tau[1,1] * W1ML)) / population[1,2])) + (rho[1,2] * ((Y1MM +
(delta[1,1] * tau[1,1] * W1MM)) / population[2,2])) + (rho[1,3] * ((Y1MH +
(delta[1,1] * tau[1,1] * W1MH)) / population[3,2]))))) - ((mew + gamma[4,1])
* W1FL);
dW1FM.dt = (VFM * (partner[2,1] * beta[1,1] * ((rho[2,1] *
((Y1ML +
(delta[1,1] * tau[1,1] * W1ML)) / population[1,2])) + (rho[2,2] * ((Y1MM +
(delta[1,1] * tau[1,1] * W1MM)) / population[2,2])) + (rho[2,3] * ((Y1MH +
(delta[1,1] * tau[1,1] * W1MH)) / population[3,2]))))) - ((mew + gamma[4,1])
* W1FM);
dW1FH.dt = (VFH * (partner[3,1] * beta[1,1] * ((rho[3,1] *
((Y1ML +
(delta[1,1] * tau[1,1] * W1ML)) / population[1,2])) + (rho[2,2] * ((Y1MM +
(delta[1,1] * tau[1,1] * W1MM)) / population[2,2])) + (rho[3,3] * ((Y1MH +
(delta[1,1] * tau[1,1] * W1MH)) / population[3,2]))))) - ((mew + gamma[4,1])
* W1FH);
dW1ML.dt = (VML * (partner[1,1] * beta[2,1] * ((rho[1,1] *
((Y1FL +
(delta[2,1] * tau[2,1] * W1FL)) / population[1,1])) + (rho[1,2] * ((Y1FM +
(delta[2,1] * tau[2,1] * W1FM)) / population[2,1])) + (rho[1,3] * ((Y1FH +
(delta[2,1] * tau[2,1] * W1FH)) / population[3,1]))))) - ((mew + gamma[3,1])
* W1MM);
dW1MM.dt = (VMM * (partner[2,1] * beta[2,1] * ((rho[2,1] *
((Y1FL +
(delta[2,1] * tau[2,1] * W1FL)) / population[1,1])) + (rho[2,2] * ((Y1FM +
(delta[2,1] * tau[2,1] * W1FM)) / population[2,1])) + (rho[2,3] * ((Y1FH +
(delta[2,1] * tau[2,1] * W1FH)) / population[3,1]))))) - ((mew + gamma[3,1])
* W1MM);
dW1MH.dt = (VMH * (partner[3,1] * beta[2,1] * ((rho[3,1] *
((Y1FL +
(delta[2,1] * tau[2,1] * W1FL)) / population[1,1])) + (rho[2,2] * ((Y1FM +
(delta[2,1] * tau[2,1] * W1FM)) / population[2,1])) + (rho[3,3] * ((Y1FH +
(delta[2,1] * tau[2,1] * W1MH)) / population[3,1]))))) - ((mew + gamma[3,1])
* W1MH);
dW2FL.dt = (VFL * (partner[1,1] * beta[1,1] * ((rho[1,1] *
((Y2ML +
(delta[1,2] * tau[1,2] * W2ML)) / population[1,2])) + (rho[1,2] * ((Y1MM +
(delta[1,2] * tau[1,2] * W2MM)) / population[2,2])) + (rho[1,3] * ((Y2MH +
(delta[1,2] * tau[1,2] * W2MH)) / population[3,2]))))) - ((mew + gamma[4,2])
* W1FL);
dW2FM.dt = (VFM * (partner[2,1] * beta[1,1] * ((rho[2,1] *
((Y2ML +
(delta[1,2] * tau[1,2] * W2ML)) / population[1,2])) + (rho[2,2] * ((Y1MM +
(delta[1,2] * tau[1,2] * W2MM)) / population[2,2])) + (rho[2,3] * ((Y2MH +
(delta[1,2] * tau[1,2] * W2MH)) / population[3,2]))))) - ((mew + gamma[4,2])
* W1FM);
dW2FH.dt = (VFH * (partner[3,1] * beta[1,1] * ((rho[3,1] *
((Y2ML +
(delta[1,2] * tau[1,2] * W2ML)) / population[1,2])) + (rho[2,2] * ((Y1MM +
(delta[1,2] * tau[1,2] * W2MM)) / population[2,2])) + (rho[3,3] * ((Y2MH +
(delta[1,2] * tau[1,2] * W2MH)) / population[3,2]))))) - ((mew + gamma[4,2])
* W1FH);
dW2ML.dt = (VML * (partner[1,1] * beta[2,1] * ((rho[1,1] *
((Y2FL +
(delta[2,2] * tau[2,2] * W2FL)) / population[1,1])) + (rho[1,2] * ((Y1FM +
(delta[2,2] * tau[2,2] * W2FM)) / population[2,1])) + (rho[1,3] * ((Y2FH +
(delta[2,2] * tau[2,2] * W2FH)) / population[3,1]))))) - ((mew + gamma[3,2])
* W1MM);
dW2MM.dt = (VMM * (partner[2,1] * beta[2,1] * ((rho[2,1] *
((Y2FL +
(delta[2,2] * tau[2,2] * W2FL)) / population[1,1])) + (rho[2,2] * ((Y1FM +
(delta[2,2] * tau[2,2] * W2FM)) / population[2,1])) + (rho[2,3] * ((Y2FH +
(delta[2,2] * tau[2,2] * W2FH)) / population[3,1]))))) - ((mew + gamma[3,2])
* W1MM);
dW2MH.dt = (VMH * (partner[3,1] * beta[2,1] * ((rho[3,1] *
((Y2FL +
(delta[2,2] * tau[2,2] * W2FL)) / population[1,1])) + (rho[2,2] * ((Y1FM +
(delta[2,2] * tau[2,2] * W2FM)) / population[2,1])) + (rho[3,3] * ((Y2FH +
(delta[2,2] * tau[2,2] * W2MH)) / population[3,1]))))) - ((mew + gamma[3,2])
* W1MH);
return(list(c(dXFL.dt, dXFM.dt, dXFH.dt, dXML.dt, dXMM.dt,
dXMH.dt,
dY1FL.dt, dY1FM.dt, dY1FH.dt, dY1ML.dt, dY1MM.dt, dY1MH.dt, dY2FL.dt,
dY2FM.dt, dY2FH.dt, dY2ML.dt, dY2MM.dt, dY2MH.dt, dZFL.dt, dZFM.dt, dZFH.dt,
dZML.dt, dZMM.dt, dZMH.dt, dVFL.dt, dVFM.dt, dVFH.dt, dVML.dt, dVMM.dt,
dVMH.dt, dW1FL.dt, dW1FM.dt, dW1FH.dt, dW1ML.dt, dW1MM.dt, dW1MH.dt,
dW2FL.dt, dW2FM.dt, dW2FH.dt, dW2ML.dt, dW2MM.dt, dW2MH.dt)));
})
}
#giving the parameters
mew = 1/15
#proportion of individuals entering
or exiting the sexually active group at a time
total = 60020
#total
population of sexually
active
phi = 0.9
#Proportion of individuals who are
successfully vaccinated
sigma = 1/10
#loss
of vaccination status
gamma = matrix(data=c(0.66, 0.66, 0.66, 0.66, 0.66, 0.66, 0.66, 0.66),
ncol=2, nrow=4)
#Duration of
infectiousness
omega = matrix(data=c(0.82, 0.15, 0.03))
#proportion
of those in each sexual activity level
population = matrix(data=c(omega[1,1]*0.5*total, omega[2,1]*0.5*total,
omega[3,1]*0.5*total, omega[1,1]*0.5*total, omega[2,1]*0.5*total,
omega[3,1]*0.5*total), ncol=2, nrow=3) #population in each sexual activity
level and gender group
partner = matrix(data=c(1.4, 3, 9))
#average number
of partners per year by risk group
beta = matrix(data=c(0.8, 0.7))
#Chance of
infection
from partner given 1 sexual encounter
rho = matrix(data=c(1, 1, 1, 1, 1, 1, 1, 1, 1), ncol=3, nrow=3)
#preference for selecting sexual partner by activity group
tau = matrix(data=c(1, 1, 1, 1), ncol=2, nrow=2)
delta = matrix(data=c(1, 1, 1, 1), ncol=2, nrow=2)
#reduction of infection from a vaccinated infected individual
XFL0 = 1850 #Initial number of
females in the low risk group
unvaccinated susceptibles
XFM0 = 340 #Initial number of
susceptible unvaccinated females medium
risk
XFH0 = 70 #Initial number of
susceptible unvaccinated females high risk
XML0 = 1850 #Initial number of
susceptible unvaccinated males low risk
XMM0 = 340 #Initial number of
susceptible unvaccinated males medium
risk
XMH0 = 70 #Initial number of
susceptible unvaccinated males high risk
Y1FL0 = 590 #Initial number of
infected unvaccinated females low risk
Y1FM0 = 100 #Initial number of
infected unvaccinated females medium
risk
Y1FH0 = 20 #Initial number of
infected unvaccinated females high risk
Y1ML0 = 590 #Initial number of
infected unvaccinated males low risk
Y1MM0 = 100 #Initial number of
infected unvaccinated males medium risk
Y1MH0 = 20 #Initial number of
infected unvaccinated males high risk
Y2FL0 = 590 #Initial number of
infected unvaccinated females low risk
Y2FM0 = 100 #Initial number of
infected unvaccinated females medium
risk
Y2FH0 = 20 #Initial number of
infected unvaccinated females high risk
Y2ML0 = 590 #Initial number of
infected unvaccinated males low risk
Y2MM0 = 100 #Initial number of
infected unvaccinated males medium risk
Y2MH0 = 20 #Initial number of
infected unvaccinated males high risk
ZFL0 = 0 #Initial number of
immune females low risk
ZFM0 = 0 #Initial number of
immune females medium risk
ZFH0 = 0 #Initial number of
immune females high risk
ZML0 = 0 #Initial number of
immune males low risk
ZMM0 = 0 #Initial number of
immune males medium risk
ZMH0 = 0 #Initial number of
immune males high risk
VFL0 = 21070 #Initial number of susceptible
vaccinated females low risk
VFM0 = 3850 #Initial number of
susceptible vaccinated females medium
risk
VFH0 = 770 #Initial number of
susceptible vaccinated females high risk
VML0 = 21070 #Initial number of susceptible
vaccinated males low risk
VMM0 = 3850 #Initial number of
susceptible vaccinated males medium risk
VMH0 = 770 #Initial number of
susceptible vaccinated males high risk
W1FL0 = 1110 #Initial number of infected
vaccinated females low risk
W1FM0 = 200 #Initial number of
infected vaccinated females medium risk
W1FH0 = 40 #Initial number of
infected vaccinated females high risk
W1ML0 = 1110 #Initial number of infected
vaccinated males low risk
W1MM0 = 200 #Initial number of
infected vaccinated males medium risk
W1MH0 = 40 #Initial number of
infected vaccinated males high risk
W2FL0 = 1110 #Initial number of infected
vaccinated females low risk
W2FM0 = 200 #Initial number of
infected vaccinated females medium risk
W2FH0 = 40 #Initial number of
infected vaccinated females high risk
W2ML0 = 1110 #Initial number of infected
vaccinated males low risk
W2MM0 = 200 #Initial number of
infected vaccinated males medium risk
W2MH0 = 40 #Initial number of
infected vaccinated males high risk
p = list(mew=mew, total=total, phi=phi, sigma=sigma, gamma=gamma,
omega=omega, population=population, partner=partner, beta=beta, rho=rho,
tau=tau, delta=delta)
y0 = c(XFL0, XFM0, XFH0, XML0, XMM0, XMH0, Y1FL0, Y1FM0, Y1FH0, Y1ML0,
Y1MM0, Y1MH0, Y2FL0, Y2FM0, Y2FH0, Y2ML0, Y2MM0, Y2MH0, ZFL0, ZFM0, ZFH0,
ZML0, ZMM0, ZMH0, VFL0, VFM0, VFH0, VML0, VMM0, VMH0, W1FL0, W1FM0, W1FH0,
W1ML0, W1MM0, W1MH0, W2FL0, W2FM0, W2FH0, W2ML0, W2MM0, W2MH0)
#Running the ode integrator
steps= 10;
t = seq(from=0, to=100, by=1);
Is anyone able to help?
--
View this message in context:
http://r.789695.n4.nabble.com/DeSolver-giving-NA-as-output-but-running-fully-tp4706497.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
