[R] DeSolver giving "NA" as output, but running fully.
John Kane
jrkrideau at inbox.com
Mon Apr 27 23:46:46 CEST 2015
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!
More information about the R-help
mailing list