thr3ads.net - R help - [R] DeSolver giving "NA" as output, but running fully. [Apr 2015]

If this information is useful, please help other people find it:
Share via:

walke554

2015-Apr-27 20:34 UTC

[R] DeSolver giving "NA" as output, but running fully.

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.

John Kane

2015-Apr-27 21:46 UTC

head link

[R] DeSolver giving "NA" as output, but running fully.

Data?

Use deput()  (see ?dput) to provide some sample data.  Also you might find this
useful http://adv-r.had.co.nz/Reproducibility.html

John Kane
Kingston ON Canada

> -----Original Message-----
> From: walke554 at umn.edu
> Sent: Mon, 27 Apr 2015 13:34:54 -0700 (PDT)
> To: r-help at r-project.org
> Subject: [R] DeSolver giving "NA" as output, but running fully.
> 
> 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 at 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.
____________________________________________________________
FREE 3D EARTH SCREENSAVER - Watch the Earth right on your desktop!

walke554

2015-Apr-29 21:27 UTC

head link

[R] DeSolver giving "NA" as output, but running fully.

All of the data you need to run the ODE integrator is there.  all the
parameters and starting populations and values are given under the:

"#giving the parameters"

section of the code.



--
View this message in context:
http://r.789695.n4.nabble.com/DeSolver-giving-NA-as-output-but-running-fully-tp4706497p4706620.html
Sent from the R help mailing list archive at Nabble.com.

Tsjerk Wassenaar

2015-Apr-30 07:01 UTC

head link

[R] DeSolver giving "NA" as output, but running fully.

Hey :)


        W2MH = y[43]; #number of infected vaccinated males high
risk> infected with
> non-vaccine strain
>
> length(y0)
[1] 42


As a sidenote, would you mind sharing the flow diagram with me, so I can
show it to my students doing a practical with DeSolve, as example of a
contemporary model?

Cheers,

Tsjerk

-- 
Tsjerk A. Wassenaar, Ph.D.

	[[alternative HTML version deleted]]

R help - Apr 2015 - DeSolver giving "NA" as output, but running fully.

[R] DeSolver giving "NA" as output, but running fully.

[R] DeSolver giving "NA" as output, but running fully.

[R] DeSolver giving "NA" as output, but running fully.

[R] DeSolver giving "NA" as output, but running fully.