From 353863992601d803574c8c4c27f0621518883b14 Mon Sep 17 00:00:00 2001 From: "Chaillet, Jack" Date: Fri, 25 Mar 2022 16:26:31 -0400 Subject: [PATCH] Added Example Time Series Code --- ResistanceTimeSeries.R | 117 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 117 insertions(+) create mode 100644 ResistanceTimeSeries.R diff --git a/ResistanceTimeSeries.R b/ResistanceTimeSeries.R new file mode 100644 index 0000000..caeb41d --- /dev/null +++ b/ResistanceTimeSeries.R @@ -0,0 +1,117 @@ +t4 = seq(from=0,to=4*365,by=1) +CQvec = c(1,1,1,1,0.54,0.54,0.54,0.54,0.5,0.5,0.5,0.5,0.18,0.18,0.18,0.18,0.01,0.01,0.01,0.01) +CQvecGhana = c(1,1,1,1,0.3,0.3,0.3,0.3,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.1,0.1,0.1,0.1) +CQvecnull = c(0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02) +CQvecmax = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1) +eqstagSU = matrix(data=0,nrow=21,ncol=15) +eqstagSD = matrix(data=0,nrow=21,ncol=15) +eqstagRU = matrix(data=0,nrow=21,ncol=15) +eqstagRD = matrix(data=0,nrow=21,ncol=15) +eqstagNU = matrix(data=0,nrow=21,ncol=15) +eqstagND = matrix(data=0,nrow=21,ncol=15) +eqstagXU = matrix(data=0,nrow=21,ncol=15) +eqstagXD = matrix(data=0,nrow=21,ncol=15) +eqstagAU = matrix(data=0,nrow=21,ncol=15) +eqstagAD = matrix(data=0,nrow=21,ncol=15) +eqosc = matrix(data=0,nrow=21,ncol=15) +eqstagf = matrix(data=0,nrow=21,ncol=15) +for (k in 1:15){ + wr=1-0.6 + bavgnull = 10^bavgvaluesverynarrow[k] + eqstagSU[1,k] = out3eqSU[k] + #eqstagSU[1,k] = (1-out3eq2[k])*1000 + eqstagSD[1,k] = out3eqSD[k] + #eqstagSD[1,k] = (1-out3eq2[k])*1000 + eqstagRU[1,k] = out3eqRU[k] + #eqstagRU[1,k] = out3eq2[k]*1000 + eqstagRD[1,k] = out3eqRD[k] + #eqstagRD[1,k] = out3eq2[k]*1000 + eqstagf[1,k] = (eqstagRU[1,k]+eqstagRD[1,k])/(eqstagSU[1,k]+eqstagSD[1,k]+eqstagRU[1,k]+eqstagRD[1,k]) + eqstagNU[1,k] = out3eqNU[k] + #eqstagNU[1,k] = 333 + eqstagND[1,k] = out3eqND[k] + #eqstagND[1,k] = 333 + eqstagXU[1,k] = out3eqXU[k] + #eqstagXU[1,k] = 333 + eqstagXD[1,k] = out3eqXD[k] + #eqstagXD[1,k] = 333 + eqstagAU[1,k] = out3eqAU[k] + #eqstagAU[1,k] = 333 + eqstagAD[1,k] = out3eqAD[k] + #eqstagAD[1,k] = 333 + eqosc[1,k] = out3osc[k] + #eqosc[1,k] = bavgnull+(bavgnull*L) + for (i in 1:20){ + CQ=CQvec[i] + #N0 = c(eqstagSU[i,k],eqstagSD[i,k],eqstagRU[i,k],eqstagRD[i,k],333,333,333,333,333,333,(bavgnull)+(bavgnull*L)) + N0 = c(eqstagSU[i,k],eqstagSD[i,k],eqstagRU[i,k],eqstagRD[i,k],eqstagNU[i,k],eqstagND[i,k],eqstagXU[i,k],eqstagXD[i,k],eqstagAU[i,k],eqstagAD[i,k],eqosc[i,k]) + p2 = list(ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att,bavgnull=bavgnull,CQ=CQ) + out8 = ode(y=N0,times=t4,func=SEAI,parms=p2) + #eqstag[i+1] = (out8[,4]+out8[,5])/(out8[,2]+out8[,3]+out8[,4]+out8[,5]) + eqstagSU[i+1,k] = out8[1*365,2] + eqstagSD[i+1,k] = out8[1*365,3] + eqstagRU[i+1,k] = out8[1*365,4] + eqstagRD[i+1,k] = out8[1*365,5] + eqstagNU[i+1,k] = out8[1*365,6] + eqstagND[i+1,k] = out8[1*365,7] + eqstagXU[i+1,k] = out8[1*365,8] + eqstagXD[i+1,k] = out8[1*365,9] + eqstagAU[i+1,k] = out8[1*365,10] + eqstagAD[i+1,k] = out8[1*365,11] + eqosc[i+1,k] = out8[1*365,12] + eqstagf[i+1,k] = (eqstagRU[i+1,k]+eqstagRD[i+1,k])/(eqstagSU[i+1,k]+eqstagSD[i+1,k]+eqstagRU[i+1,k]+eqstagRD[i+1,k]) + } +} +plot(eqstagf[,1],ylim=c(0,1),type='l',xlim=c(0,21)) +lines(eqstagf[,2]) +lines(eqstagf[,3]) +lines(eqstagf[,4]) +lines(eqstagf[,5]) +lines(eqstagf[,6]) +lines(eqstagf[,7]) +lines(eqstagf[,8]) +lines(eqstagf[,9]) +lines(eqstagf[,10]) +lines(eqstagf[,11]) +lines(eqstagf[,12]) +lines(eqstagf[,13]) +lines(eqstagf[,14]) +lines(eqstagf[,15]) +abline(h=1,col="BLUE") +abline(h=mean((WWARN_ACT_Partner_Drug_Mol_Surveyor_Data$present)/(WWARN_ACT_Partner_Drug_Mol_Surveyor_Data$tested)),col="RED") +abline(h=mean((`WWARN_ACT_Partner_Drug_Mol_Surveyor_Data(1)`$present)/(`WWARN_ACT_Partner_Drug_Mol_Surveyor_Data(1)`$tested)),col="GREEN") +plot(WWARN_ACT_Partner_Drug_Mol_Surveyor_Data$study.start-1998,(WWARN_ACT_Partner_Drug_Mol_Surveyor_Data$present)/(WWARN_ACT_Partner_Drug_Mol_Surveyor_Data$tested)) +plot(`WWARN_ACT_Partner_Drug_Mol_Surveyor_Data(1)`$study.start-2005,(`WWARN_ACT_Partner_Drug_Mol_Surveyor_Data(1)`$present)/(`WWARN_ACT_Partner_Drug_Mol_Surveyor_Data(1)`$tested)) + +plot(analysis_PointsGhana$year,analysis_PointsGhana$value) + +lines(eqstagf[,16]) +lines(eqstagf[,17]) +lines(eqstagf[,18]) +lines(eqstagf[,19]) +lines(eqstagf[,20])#,col="RED") +lines(eqstagf[,21]) +lines(eqstagf[,22])#,col="GREEN") +lines(eqstagf[,23]) + + +times = matrix(data=0,nrow=21,ncol=23) +for (i in 1:21){ + for(j in 1:23){ + times[i,j] = i + } +} +print(times) +timesvec = c(times) +print(timesvec) + + +trajectory = eqstagf +trajectoryvec = c(trajectory) +trajectory2 = cbind(timesvec,vec7t2,trajectoryvec) +trajectories.df <- as.data.frame(trajectory2)#,row.names = columns, col.names=rows) +rownames(trajectories.df) +times=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20) + + +ggplot(data=trajectories.df,mapping=aes(x=timesvec,y=trajectoryvec))+geom_line(aes(color=as.factor(vec7t2)))#+xlim(2,21)