ageSPVL_annals.Rmd

---
title: "Overview record of the AgeAndSPVL project"
output:
  html_document:
    df_print: paged
    theme: cosmo
    css: tablestyle.css
  pdf_document: default
  word_document: default
---

## Introduction

This document is initiated on August 4, 2018. However, the project itself was begun in late January and early February, 2018, but then progress was halted again until late July. Steve set up a Github repository back in February, but did not commit most of the files from that period until just now. The early period was also highly exploratory, with no documentation.  This is an attempt to document and standardize the work so far and then document all subsequent work moving forward.

## Original questions:

Steve asked: It has been observed that SPVL increases with age at seroconversion, although it is not well understood why. Given the results of the role paper, when EI men had higher mSPVL than ER/RV men because of the higher selective pressure caused by the narrowed transmission bottleneck, he hypothesized that perhaps something similar is going on with age.  That is, older people have less opportunity for acquisition on average (fewer acts and/or partners), so that those who did get infected would do so with higher SPVL. Indeed, the fact that people are generally assortative by age seemed like it could magnify this effect, in contrast to the role effect where the alternating chains of infection (I->R->I->R) kept the two from diverging too far.

When Steve brought up looking at age and SPVL, John had a different question: do populations in which the risk is heavily concentrated in young people have overall higher mSPVL than populations with the same overall amount of risk but spread out more over the lifecourse?

## Folder structure

The top-level `AgeAndSPVL` directory contains the .R files that each run a different scenario, and the corresponding .rda and .pef files that contain output. At this point, the names of these files are:

`ageSPVL_mXX.R`
`agePSVL_mXX.rda`
`agePSVL_mXX.pdf`

For some scenarios there are also files created that print the viral load trajectories of all infectd agents as part of understanding the underlying dynamics.

The directory also contains a set of files beginning with the name `ageSPVL_annals`, including this .Rmd file and its outputs.

## Initial models

Steve began by exploring some small, simple scenarios with just a single replicate in order to get some intution.  Initially (i.e. back in January 2018) there were 10 scenarios (`ageSPVL_mXX` where XX = 01 to 10). originally these were compiled and analyzed in the file ageSPVL_explore_meta.R, but that code is now subsumed below.

Parameters looked at in these intial runs were:

Run | 1	| 2	| 3	| 4	| 5	| 6	| 7	| 8	| 9	| 10
---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- 
min age	| 18 | 18 | 18 | 18 | 18 | 18 | 18 | 18 | 18 | 18
maxage | 55 | 55 | 55 | 55 | 55 | 55 | 55 | 55 | 55 | 55
mean_sqrtage_diff | 0.3 | 0.3 | 0.3 | 0.3 | 0.3 | 0.3 | 1 | 1.5 | 2 | 0.6
mean_sex_acts_per_day | 0.2 | --- | --- | --- | --- | --- | --- | --- | --- | ---
prob_sex_by_age | F | T | T | T | T | T | T | T | T | T
prob_sex_age_19 | --- | 0.4 | 0.3 | 0.2 | 0.4 | 0.4 | 0.3 | 0.3 | 0.3 | 0.3
max_age_sex | --- | 55 | 55 | 55 | 35 | 55 | 55 | 55 | 55 | 35
relation_dur | 50 | 50 | 50 | 50 | 50 | 200 | 50 | 50 | 50 | 50

Compiling runs:
```{R m01-m10}
nruns <- 10
popatts <- popsumm <- list()

agecoef <- iSPVL <- ageinf <- prev <- numinc <- agematch <- rep(NA, nruns)

for (i in 1:nruns) {
  if(i<10) filler <- "0" else filler <- ""
  load(paste("experiments/msm/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  popatts[[i]] <- obj$pop[[1]]
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm[[1]]
  agecoef[i] <- lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])$coef[[2]]
  iSPVL[i] <- mean(log10(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0]),na.rm=TRUE)
  ageinf[i] <- mean(popatts[[i]]$age_infection,na.rm=TRUE)
  prev[i] <- tail(popsumm[[i]]$prevalence,1)
  numinc[i] <- sum(popatts[[i]]$Time_Inf>0, na.rm=TRUE)
  agematch[i] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
  rm(obj)
}
```


Analysis:

```{R}
plot(numinc, prev, type='b')
```

Incidence and prevalence are nearly perfectly correlated (no surprise, just a good check).

```{R}
plot(agematch)
```

Mean square root age difference of 1 is about what is expected by chance in this model.

```{R}
plot(prev, iSPVL, type='b')
```

Prevalence generally predicts average SPVL.  However, runs 8, 9 and 10 generally lie above the relational line. Runs 8 and 9 are actually the cases that had disassortative mixing by age (bigger difference than expected by chance).  And 10 concentrates the sex at younger ages and has moderate assortativity by age. So now we're getting somewhere.

```{R}
plot(agecoef); abline(h=0)
```

The only cases where the coefficient on SPVL~age is positive are the ones where there is disassortative mixing by age.  The others are the reverse.  This is very very strange and requires much probing.

Some thoughts Steve has at this point: perhaps these runs are far from equilibrium.  Does getting infected younger mean getting infected later in the simulation on average, and there is a secular trend in SPVL increasing overall?

This also raises the question of whether this paper should also think about trends in SPVL over the course of the epidemic finally.

Next step: run regression with date of infection and age of infection to see what happens.

## Aug 6, 2018

Ran the regression with date infected as well:

```{R}
agecoef2 <- rep(NA, nruns)

for (i in 1:nruns) {
  agecoef2[i] <- lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
    popatts[[i]]$Time_Inf[popatts[[i]]$Time_Inf>0]+
    popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])$coef[[3]]
}
plot(agecoef2)

```

Also, curious to confirm the ages at which infection is occuring:

```{R, out.width="20%"}

meanageinf <- rep(NA,10)
for (i in 1:nruns) {
    hist(popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0],main="",
        breaks=c(15,20,25,30,35,40,45,50,55,60))
    meanageinf[i] <- mean(popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0], na.rm=TRUE)
}

```

## Aug 8, 2018

Removed agePSVL-explore_meta.R because I've realized that it is redundant with this document -- the code resides here as a living document.

The distribution of ages at infection in runs 1 through 10 show that only runs 8 and 9 have a sizeable number of people getting infected at older ages. And they're the two with SPVL increasing with age.  They're also the two that have age-discordant mixing.  But maybe it's not the age-discordant mixing per se that creates the positive age/SPVL effect -- maybe it's just having some  older infections.  So some new runs:

Run | 11 | 12 | 13
---- | ---- | ---- | ----
min age	| 18 | 18 | 18
maxage | 55 | 55 | 55
mean_sqrtage_diff | 0.6 | 0.6 | 0.6
mean_sex_acts_per_day | --- | --- | ---
prob_sex_by_age | T | T | T
prob_sex_age_19 | 0.2 | 0.2 | 0.2
max_age_sex | 55 | 55 | 55
relation_dur | 50 | 100 | 200

```{R m11-m13}
for (i in 11:13) {
  if(i<10) filler <- "0" else filler <- ""
  load(paste("experiments/msm/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  popatts[[i]] <- obj$pop[[1]]
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm[[1]]
  agecoef[i] <- lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])$coef[[2]]
  iSPVL[i] <- mean(log10(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0]),na.rm=TRUE)
  ageinf[i] <- mean(popatts[[i]]$age_infection,na.rm=TRUE)
  prev[i] <- tail(popsumm[[i]]$prevalence,1)
  numinc[i] <- sum(popatts[[i]]$Time_Inf>0, na.rm=TRUE)
  agematch[i] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
  meanageinf[i] <- mean(popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0], na.rm=TRUE)
  rm(obj)
}

plot(numinc, prev, type='b')
plot(prev, iSPVL, type='b')
plot(agecoef); abline(h=0)
plot(agematch)
```

So run 12 is maybe above?  Which is odd - it's not monotonoic (or it's just stochastic).

Let's look at mean age infected as a summary stat on the age dist:

```{R}
plot(meanageinf, agecoef)
```

Let's redo 11-13 as 14-16 with higher overall incidence so that we can have more data:


Run | 14 | 15 | 16
---- | ---- | ---- | ----
min age	| 18 | 18 | 18
maxage | 55 | 55 | 55
mean_sqrtage_diff | 0.6 | 0.6 | 0.6
mean_sex_acts_per_day | --- | --- | ---
prob_sex_by_age | T | T | T
prob_sex_age_19 | 0.3 | 0.3 | 0.3
max_age_sex | 55 | 55 | 55
relation_dur | 50 | 100 | 200

```{R m14-m16}
for (i in 14:16) {
  if(i<10) filler <- "0" else filler <- ""
  load(paste("experiments/msm/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  popatts[[i]] <- obj$pop[[1]]
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm[[1]]
  agecoef[i] <- lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])$coef[[2]]
  iSPVL[i] <- mean(log10(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0]),na.rm=TRUE)
  ageinf[i] <- mean(popatts[[i]]$age_infection,na.rm=TRUE)
  prev[i] <- tail(popsumm[[i]]$prevalence,1)
  numinc[i] <- sum(popatts[[i]]$Time_Inf>0, na.rm=TRUE)
  agematch[i] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
  meanageinf[i] <- mean(popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0], na.rm=TRUE)
  rm(obj)
}

plot(numinc, prev, type='b')
plot(prev, iSPVL, type='b')
plot(agecoef); abline(h=0)
plot(agematch)
plot(meanageinf)
plot(meanageinf, agecoef, type='b')
plot(meanageinf, iSPVL)
```

The evidence for John's hypothesis is beginning to accumulate.
For Steve's, not so much.

Steve is pondering more. Is older age really like being insertive?  That is, insertive guys have as many partners and exposures as receptive guys do (actually, more exposures given higher prevalence in their partner pool); but they have lower probability of acquisition per act.  We're modeling older people as having fewer acts than others.  Is that really an analogous mechanism? 

Let's dive deeper into two runs, 12 and 15:

```{R}
p12 <- popatts[[12]]
p15 <- popatts[[15]]
plot(p12$age_infection, p12$Donors_age)
plot(p15$age_infection, p15$Donors_age)
plot(p12$Donors_age, p12$Donors_LogSetPoint)
plot(p15$Donors_age, p15$Donors_LogSetPoint)
```

Aha. Steve has a sudden realization, of something that he thinks Sarah actually hypothesized back in Feb/Mar, but which he forgot about until now.

Because there is no treatment, people with high SPVL die quickly.  With assortative age mixing, older folks mostly get infected by other older folks, but most of the folks with high SPVL have died before they reach older age.  So there are two different effects operating here in different directions.

The trick will be to add treatment into the model, in ways that make it so that when treatment fails they go back up to the SPVL they would have in the absence of treatment.

So, things Steve needs to figure out:
(1) how do the VL dynamics work in the presence of treatment
(2) also - how does the prob_sex_by_age / prob_sex_age_19  / max_age_sex code work exactly?  Why can't max_age_sex be > 55?

Steve is also thinking about an additional analysis: if one has heterosexual asymmetric age mixing (a la absdiffby) does that lead to higher mSPVL for women, even after controlling for age at infection?  This could be an interesting analysis, but also requires understanding the lit on sex differences in SPVL a bit more - does it appear in setting with tx, without tx, or both?

OK, as an experiment, trying a version in which the risk is concentrated exclusively in a 10-year period, but still where risk declines over that period. These folks shouldn't see much die-off among partners.  Let's see....

Run | 17 | 18 | 19
---- | ---- | ---- | ----
min age	| 18 | 18 | 18
maxage | 55 | 55 | 55
mean_sqrtage_diff | 0.6 | 0.6 | 0.6
mean_sex_acts_per_day | --- | --- | ---
prob_sex_by_age | T | T | T
prob_sex_age_19 | 1 | 1 | 1
max_age_sex | 29 | 29 | 29
relation_dur | 50 | 100 | 200


```{R m17-m19}
for (i in 17:19) {
  if(i<10) filler <- "0" else filler <- ""
  load(paste("experiments/msm/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  popatts[[i]] <- obj$pop[[1]]
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm[[1]]
  agecoef[i] <- lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])$coef[[2]]
  iSPVL[i] <- mean(log10(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0]),na.rm=TRUE)
  ageinf[i] <- mean(popatts[[i]]$age_infection,na.rm=TRUE)
  prev[i] <- tail(popsumm[[i]]$prevalence,1)
  numinc[i] <- sum(popatts[[i]]$Time_Inf>0, na.rm=TRUE)
  agematch[i] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
  meanageinf[i] <- mean(popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0], na.rm=TRUE)
  rm(obj)
}

plot(numinc, prev, type='b')
plot(prev, iSPVL, type='b')
plot(agecoef); abline(h=0)
plot(agematch)
plot(meanageinf)
plot(meanageinf, agecoef, type='b')
plot(meanageinf, iSPVL)
```


```{R}
i <- 17
summary(lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0]))
```

```{R}
for (i in 17:19) {
    hist(popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0],main="",
        breaks=17:65)
}
```

```{R}

plot(popatts[[17]]$age_infection[popatts[[i]]$Time_Inf>0],
     popatts[[17]]$LogSetPoint[popatts[[i]]$Time_Inf>0]
     
       )
```

```{R}
i <- 17
summary(lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0 & popatts[[i]]$age_infection<30],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0 & popatts[[i]]$age_infection<30]))
```

## Aug 9-Sep 20 2018

I learned from John and James that they have code to model tx cessation and VL rebound. However, the arguments for these are a little unclear. I spent some time working through this all in the period Aug 9-14.  However, I didn't document it as I went, leaving me to now recreate it a month later. Lesson to be learned here!!

Run 20 was my first attempt, which did not work - once people became suppressed, they stayed there; or, if they lost supression, they immediately regained it.

From there, James and John and I shared a series of emails (subject: treatment dropout scripts/results, dates Aug 9-11). James revealed some arguments that needed to be set, and the fact that the "prob_tx_droput" argument is spelled wrong (missing an o). He also included a sample script that worked, which I have included as ageSPVL_m21_jtm2.R. This had some duplicate arguments in it (ones set twice), and some that didn't have to be set (including different stages of a treatment campaign). I then spent time trying to streamline it while still getting the behavior to work as expected (i.e. people having both long periods of suppression and then subsequent periods of non-suppression.) One of these attempts is ageSPVL_m21.R. Then James realized he needed to commit the function targeted_treatment2 to the master branch, which has the needed code to assign treatment randomly.  Runs 22-25 use this, and explore different scenarios, although they do not appear to be all that different from one another. However, I did not finish ascertaining what every parameter does and which ones are truly needed, nor did  I determine what proportion of people are actually on treatment at any point in time, what proportion go onto treatment immediately, etc.  This is all crucial for being able to understand the scenarios being modeled and relate observed outcomes to their features.  Looking at it now I see the following results, although I'm not ready to put much stock in their interpretation:

Run | 22 | 23 | 24 | 25
---- | ---- | ---- | ---- | ----
min age	| 18 | 18 | 18 | 18
maxage | 55 | 55 | 55 | 55
mean_sqrtage_diff | 0.6 | 0.6 | 0.3 | 0.3
mean_sex_acts_per_day | --- | --- | --- | ---
prob_sex_by_age | T | T | T | T
prob_sex_age_19 | 0.4 | 0.4 | 0.4 | 0.4 
max_age_sex | 55 | 55 | 55 | 55
relation_dur | 200 | 200 | 200 | 200 
tx_type | "random" | "random" | "random" | "random" 
mean_trtmnt_delay | 0 | 0 | 0 | 0
start_treatment_campaign | 1 |  1 | 1 | 1
proportion_treated | 0.1 | 0.15 | 0.15 | 0.17
testing_model | "interval" | "interval" | "interval" | "interval" 
mean_test_interval_male | 365 | 365 | 365 | 365
prob_tx_dropout | 0.1 | 0.1 | 0.1 | 0.1 


```{R m21-m25}
for (i in 21:25) {
  if(i<10) filler <- "0" else filler <- ""
  load(paste("experiments/msm/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  popatts[[i]] <- obj$pop[[1]]
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm[[1]]
  agecoef[i] <- lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])$coef[[2]]
  iSPVL[i] <- mean(log10(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0]),na.rm=TRUE)
  ageinf[i] <- mean(popatts[[i]]$age_infection,na.rm=TRUE)
  prev[i] <- tail(popsumm[[i]]$prevalence,1)
  numinc[i] <- sum(popatts[[i]]$Time_Inf>0, na.rm=TRUE)
  agematch[i] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
  meanageinf[i] <- mean(popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0], na.rm=TRUE)
  rm(obj)
}
```

```{R}
plot(numinc, prev, type='b')
```

This demonstrates that treatment is certainly working in some form, given the strong distinction between the tx and no-tx runs in this relationship.

```{R}
plot(prev, iSPVL, type='b')
plot(agecoef); abline(h=0)
plot(agematch)
plot(meanageinf)
plot(meanageinf, agecoef, type='b')
plot(meanageinf, iSPVL)
```

In the next few plots I look into whether/how the VL code works:

```{R}
plot(log10(ageSPVL_m22$vl_list[[1]][[7200]][,'vl']))
```

```{R}
plot(popatts[[23]]$age_infection, popatts[[23]]$Donors_age)
plot(popatts[[24]]$age_infection, popatts[[24]]$Donors_age)
```

```{R}
table(popatts[[24]]$Status)
```

Status | Meaning
--- | ---
1 | alive and HIV+
 2 | alive and HIV-
-1 | died of background mortality
-1.5 | aged out
-2 | died of AIDS

## October 1, 2018

I am searching through the code to see what each relevant parameter actually does:

Parameter | Use | Value w/ notes
--- | ----
start_treat_before_big_campaign | Start time of gradual ramp-up prior to the start of the big TasP campaign (this part random) | 
start_treatment_campaign | start time of tx campaign (this part depends on tx_type)
tx_type | how treatment is allocated during the tx campaig. random = anyone not_curr_tx, which means not on tx, diagnosed, and beyond mean_trtmnt_delay since diagnosis
tx_limit | flag with values "absolute_num" or "percentage"
mean_trtmnt_delay | Delay between diagnosis and availability for tx 
proportion_treated | See notes below
prob_eligible_ART | 
tx_schedule_props | % of people who always (F), sometimes (V), or never (N) take therapy
prob_tx_droput | NO "O"! Prob. of "V"-type agent discontinuing therapy over the course of a year
prop_tx_before | 
yearly_incr_tx | annual inc. in # of people being treated after the TasP campaign
prob_care | % that could get treated given "all out" campaign
prob_eligible_ART | 
vl_full_supp | 
proportion_treated_begin | Treated before ramp-up 


## October 19

Work on the above table made me realize that the code is written in a funny way, such that:

If tx_limit == "absolute_num" then on the time steps prior to the start of the tx campaign, max_num_treated is set as proportion_treated*number_infected.  Then starting with the first day of the campaign, that number becomes frozen.  Future increases in the number who can be treated is determined by yearly_incr_tx.

If tx_limit == "percentage", then the tx limit is proportion_treated*total_alive (not total_infected).

Because I had had start_treatment_campaign == 1, there was never a chance for the code to set the maximum number for treatment.  So for run 26 I changed start_treatment_campaign to 2, and proportion_treated to 0.5.  Initial run looked very promising:

```{R m26}
for (i in 26:26) {
  if(i<10) filler <- "0" else filler <- ""
  load(paste("experiments/msm/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  popatts[[i]] <- obj$pop[[1]]
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm[[1]]
  agecoef[i] <- lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])$coef[[2]]
  iSPVL[i] <- mean(log10(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0]),na.rm=TRUE)
  ageinf[i] <- mean(popatts[[i]]$age_infection,na.rm=TRUE)
  prev[i] <- tail(popsumm[[i]]$prevalence,1)
  numinc[i] <- sum(popatts[[i]]$Time_Inf>0, na.rm=TRUE)
  agematch[i] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
  meanageinf[i] <- mean(popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0], na.rm=TRUE)
  rm(obj)
}
```

```{r}
plot(agecoef); abline(h=0)
```

However, a review of the vl_traj pdf revelaed very odd things.  Done for the day.

## Oct 23, 2018

I discovered that I had tx_type set to "percentage" instead instead of "absolute_num". Changing that and upping the percentage treated to 0.5 led to behavior in the vl_traj that I would expect.  In doing this I edited run 26 instead of creating a new run, so the run 26 shown up above no longer matches what it used to be up there.

```{r}
plot(agecoef); abline(h=0)
plot(popatts[[26]]$age_infection, popatts[[26]]$Donors_age)
```

OK it's late at night and instead of finish documenting I decided to do a run 27 where I made the relationships longer, the age mixing tighter, and the probability of tx dropout lower.  Bad scientist :-)

```{r m27}
for (i in 27:27) {
  if(i<10) filler <- "0" else filler <- ""
  load(paste("experiments/msm/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  popatts[[i]] <- obj$pop[[1]]
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm[[1]]
  agecoef[i] <- lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])$coef[[2]]
  iSPVL[i] <- mean(log10(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0]),na.rm=TRUE)
  ageinf[i] <- mean(popatts[[i]]$age_infection,na.rm=TRUE)
  prev[i] <- tail(popsumm[[i]]$prevalence,1)
  numinc[i] <- sum(popatts[[i]]$Time_Inf>0, na.rm=TRUE)
  agematch[i] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
  meanageinf[i] <- mean(popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0], na.rm=TRUE)
  rm(obj)
}
```

```{r}
plot(agecoef); abline(h=0)
plot(popatts[[27]]$age_infection, popatts[[27]]$Donors_age)
i <- 27; summary(lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0]))
```
```{r}
plot(agecoef); abline(h=0)
plot(popatts[[27]]$age_infection, popatts[[27]]$Donors_age)
```

## Oct 26

Run 27 is a dud. Of course I changed three things in doing it, so it's hard to know what any of the effects are.  So I just created the runs 28-30 which change each of one of those three things back, to see if there is any clear signal across any of them. Set them running and then heading out.


```{r m28-m30}
for (i in 28:30) {
  if(i<10) filler <- "0" else filler <- ""
  load(paste("experiments/msm/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  popatts[[i]] <- obj$pop[[1]]
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm[[1]]
  agecoef[i] <- lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])$coef[[2]]
  iSPVL[i] <- mean(log10(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0]),na.rm=TRUE)
  ageinf[i] <- mean(popatts[[i]]$age_infection,na.rm=TRUE)
  prev[i] <- tail(popsumm[[i]]$prevalence,1)
  numinc[i] <- sum(popatts[[i]]$Time_Inf>0, na.rm=TRUE)
  agematch[i] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
  meanageinf[i] <- mean(popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0], na.rm=TRUE)
  rm(obj)
}
```

```{r}
plot(agecoef); abline(h=0)
plot(popatts[[28]]$age_infection, popatts[[28]]$Donors_age)
plot(popatts[[29]]$age_infection, popatts[[29]]$Donors_age)
plot(popatts[[30]]$age_infection, popatts[[30]]$Donors_age)

```

```{r}
for(i in 28:30) print(summary(lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])))
```

### Oct 30, 2018

What isn't obvious from the above results is that I ran the run 29 a few times, because its first iteration showed a high positive correlation between age and SPVL, and I wanted to see if this was robust.  My later runs reversed that.  So it seems as if the time may have come to amplify the population sizes and/or run times and/or number of runs on these to really clarify the magnitude of the effects.  Until now we wanted a "quick look" but it's now clear that, with things working we need more to know things for sure.

Note to self: still want to go and do final confirmation on how the treatment is working, and whether the expected number of people really are on treatment at any point in time.  We may also wish to do a check of how at any point in time what the viral load of those alive is by their current age (rather than their age at acquisition) because this will help to determine the magnitude of the different effects.

Note to self: Sarah wrote with a clever idea - that underlying heterogeneity in risk acquisition (whether biological or something else) which varies across person but not necessarily over time within person might also create this kind of effect. Because those who are infected later would disproportinately be those with lower per-act acquisition risk.  Should explore this too.

So, beginning now to set running three more robust runs. Starting with 28->31.  Time goes 20->50 years, nsims 1->5. (Increasing pop size leads to an error that seems to have something to do with MaxDyadTypes)

```{r m31}
for (i in 31:31) {
  if(i<10) filler <- "0" else filler <- ""
  load(paste("experiments/msm/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  popatts[[i]] <- obj$pop[[1]]
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm[[1]]
  agecoef[i] <- lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])$coef[[2]]
  iSPVL[i] <- mean(log10(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0]),na.rm=TRUE)
  ageinf[i] <- mean(popatts[[i]]$age_infection,na.rm=TRUE)
  prev[i] <- tail(popsumm[[i]]$prevalence,1)
  numinc[i] <- sum(popatts[[i]]$Time_Inf>0, na.rm=TRUE)
  agematch[i] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
  meanageinf[i] <- mean(popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0], na.rm=TRUE)
  rm(obj)
}
```

```{r}
plot(agecoef); abline(h=0)
plot(popatts[[31]]$age_infection, popatts[[31]]$Donors_age)

```

```{r}
for(i in 28:31) print(summary(lm(log(popatts[[i]]$SetPoint[popatts[[i]]$Time_Inf>0],10)~
               popatts[[i]]$age_infection[popatts[[i]]$Time_Inf>0])))
```

Standard error did go down by about half.  But it appears only one sim was done.  Perhaps this is because I didn't set ncores? Will try that now, as run 32, and bump number of runs up to 10 whole I'm at it.

Here begins a change to all data storage since there is now >1 simulation per scenario. Things that are already lists (popatts and popsumm) can remain so, but things that are vectors must be remade as lists.


```{r}
agecoef.list <- iSPVL.list <- ageinf.list <- prev.list <- 
  numinc.list <- agematch.list <- meanageinf.list <- list()
```

```{r m32}
i <- 32

if(i<10) filler <- "0" else filler <- ""
load(paste("../AgeAndSPVL_oversize/ageSPVL_m",filler,i,".rda",sep=""))
obj <- get(paste("ageSPVL_m",filler,i,sep=""))

popatts[[i]] <- obj$pop
popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm

agecoef.list[[i]] <- iSPVL.list[[i]] <- ageinf.list[[i]] <- prev.list[[i]] <- 
  numinc.list[[i]] <- agematch.list[[i]] <- meanageinf.list[[i]] <- vector()

for (j in 1:length(popatts[[i]])) {
  agecoef.list[[i]][j] <- lm(log(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0],10)~
               popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0])$coef[[2]]
  iSPVL.list[[i]][j] <- mean(log10(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0]),na.rm=TRUE)
  ageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection,na.rm=TRUE)
  prev.list[[i]][j] <- tail(popsumm[[i]][[j]]$prevalence,1)
  numinc.list[[i]][j] <- sum(popatts[[i]][[j]]$Time_Inf>0, na.rm=TRUE)
  agematch.list[[i]][j] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
  meanageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0], na.rm=TRUE)
}
rm(obj)
```

```{r}
boxplot(agecoef.list)
boxplot(iSPVL.list)
boxplot(ageinf.list)
boxplot(prev.list)
boxplot(numinc.list)
boxplot(agematch.list)
boxplot(meanageinf.list)

```

PS I seem to have the syntax for this down now. So the mapping is 28:30 -> 32:34, but with nsims = ncores = 10, and duration = 50 years.

```{r m33-m34}
for (i in 33:34) {

  if(i<10) filler <- "0" else filler <- ""
  load(paste("../AgeAndSPVL_oversize/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  
  popatts[[i]] <- obj$pop
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm
  
  agecoef.list[[i]] <- iSPVL.list[[i]] <- ageinf.list[[i]] <- prev.list[[i]] <- 
    numinc.list[[i]] <- agematch.list[[i]] <- meanageinf.list[[i]] <- vector()
  
  for (j in 1:length(popatts[[i]])) {
    agecoef.list[[i]][j] <- lm(log(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0],10)~
                 popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0])$coef[[2]]
    iSPVL.list[[i]][j] <- mean(log10(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0]),na.rm=TRUE)
    ageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection,na.rm=TRUE)
    prev.list[[i]][j] <- tail(popsumm[[i]][[j]]$prevalence,1)
    numinc.list[[i]][j] <- sum(popatts[[i]][[j]]$Time_Inf>0, na.rm=TRUE)
    agematch.list[[i]][j] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
    meanageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0], na.rm=TRUE)
  }
  rm(obj)
}
```

```{r}
boxplot(agecoef.list); abline(h=0)
boxplot(iSPVL.list)
boxplot(ageinf.list)
boxplot(prev.list)
boxplot(numinc.list)
boxplot(agematch.list)
boxplot(meanageinf.list)
```

### Nov 29, 2018

Picking up after a long time. Just went back and read through the whole long history. Notes for this moment:
- Still need to check whether the treatment works as expected! With parameters trimmed down.
- Make table for runs 32-34.
- Very interesting how strongly negative the agematch coef is for these (-3 for run 32, -10 for 33-34). Since we've seen (and it makes sense that) less age-assortative mixing can increase the possiblity for a positive result, we may have room here to lessen the age mixing.
- Am I reading that right, that prev is 4% for run 32, 2% for run 33, and 0% for run 34?


Run | 32 | 33 | 34 
---- | ---- | ---- | ---- 
min age	| 18 | 18 | 18 
maxage | 55 | 55 | 55 
mean_sqrtage_diff | *0.3* | 0.1 | 0.1 
mean_sex_acts_per_day | --- | --- | --- 
prob_sex_by_age | T | T | T 
prob_sex_age_19 | 0.4 | 0.4 | 0.4 
max_age_sex | 55 | 55 | 55 
relation_dur | 500 | *200* | 500 
tx_type | "random" | "random" | "random"
mean_trtmnt_delay | 0 | 0 | 0 
start_treatment_campaign | 1 |  1 | 1 
proportion_treated | 0.5 | 0.5 | 0.5 
testing_model | "interval" | "interval" | "interval"
mean_test_interval_male | 365 | 365 | 365
prob_tx_dropout | 0.05 | 0.05 | *0.1*

OK, now I am going to start four more runs that are:

Run | 35 | 36 | 37 | 38
---- | ---- | ---- | ---- | ---- 
min age	| 18 | 18 | 18 | 18
maxage | 55 | 55 | 55 | 55 
mean_sqrtage_diff | *0.5* | *0.5* | *0.7* | *0.7* 
mean_sex_acts_per_day | --- | --- | --- | ---
prob_sex_by_age | T | T | T | T
prob_sex_age_19 | 0.4 | 0.4 | 0.4 | 0.4
max_age_sex | 55 | 55 | 55 | 55
relation_dur | *500* | *200* | *500* | *200* 
tx_type | "random" | "random" | "random" | "random"
mean_trtmnt_delay | 0 | 0 | 0 | 0
start_treatment_campaign | 1 | 1 | 1 | 1
proportion_treated | 0.5 | 0.5 | 0.5| 0.5
testing_model | "interval" | "interval" | "interval" | "interval"
mean_test_interval_male | 365 | 365 | 365 | 365
prob_tx_dropout | 0.05 | 0.05 | 0.05 | 0.05

```{r m35-m38}
for (i in 35:38) {

  if(i<10) filler <- "0" else filler <- ""
  load(paste("../AgeAndSPVL_oversize/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  
  popatts[[i]] <- obj$pop
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm
  
  agecoef.list[[i]] <- iSPVL.list[[i]] <- ageinf.list[[i]] <- prev.list[[i]] <- 
    numinc.list[[i]] <- agematch.list[[i]] <- meanageinf.list[[i]] <- vector()
  
  for (j in 1:length(popatts[[i]])) {
    agecoef.list[[i]][j] <- lm(log(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0],10)~
                 popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0])$coef[[2]]
    iSPVL.list[[i]][j] <- mean(log10(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0]),na.rm=TRUE)
    ageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection,na.rm=TRUE)
    prev.list[[i]][j] <- tail(popsumm[[i]][[j]]$prevalence,1)
    numinc.list[[i]][j] <- sum(popatts[[i]][[j]]$Time_Inf>0, na.rm=TRUE)
    agematch.list[[i]][j] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
    meanageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0], na.rm=TRUE)
  }
  rm(obj)
}
```

```{r}
boxplot(agecoef.list); abline(h=0)
boxplot(iSPVL.list)
boxplot(ageinf.list)
boxplot(prev.list)
boxplot(numinc.list)
boxplot(agematch.list)
boxplot(meanageinf.list)
```

So run 37 isn't quite there, but it's getting close,  And the pattern may be: less assortative mixing, higher agecoef; and longer relationships, higher agecoef.  These make sense.  So let me try some more runs just to see.

Run | 39 | 40 | 41 | 42
---- | ---- | ---- | ---- | ---- 
min age	| 18 | 18 | 18 | 18
maxage | 55 | 55 | 55 | 55 
mean_sqrtage_diff | *0.7* | *0.9* | *0.9* | *1.2* 
mean_sex_acts_per_day | --- | --- | --- | ---
prob_sex_by_age | T | T | T | T
prob_sex_age_19 | 0.4 | 0.4 | 0.4 | 0.4
max_age_sex | 55 | 55 | 55 | 55
relation_dur | *1000* | *500* | *1000* | *1000* 
tx_type | "random" | "random" | "random" | "random"
mean_trtmnt_delay | 0 | 0 | 0 | 0
start_treatment_campaign | 1 | 1 | 1 | 1
proportion_treated | 0.5 | 0.5 | 0.5| 0.5
testing_model | "interval" | "interval" | "interval" | "interval"
mean_test_interval_male | 365 | 365 | 365 | 365
prob_tx_dropout | 0.05 | 0.05 | 0.05 | 0.05

```{r m39-m42}
for (i in 39:42) {

  if(i<10) filler <- "0" else filler <- ""
  load(paste("../AgeAndSPVL_oversize/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  
  popatts[[i]] <- obj$pop
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm
  
  agecoef.list[[i]] <- iSPVL.list[[i]] <- ageinf.list[[i]] <- prev.list[[i]] <- 
    numinc.list[[i]] <- agematch.list[[i]] <- meanageinf.list[[i]] <- vector()
  
  for (j in 1:length(popatts[[i]])) {
    agecoef.list[[i]][j] <- lm(log(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0],10)~
                 popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0])$coef[[2]]
    iSPVL.list[[i]][j] <- mean(log10(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0]),na.rm=TRUE)
    ageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection,na.rm=TRUE)
    prev.list[[i]][j] <- tail(popsumm[[i]][[j]]$prevalence,1)
    numinc.list[[i]][j] <- sum(popatts[[i]][[j]]$Time_Inf>0, na.rm=TRUE)
    agematch.list[[i]][j] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
    meanageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0], na.rm=TRUE)
  }
  rm(obj)
}
```

```{r}
boxplot(agecoef.list); abline(h=0)
boxplot(iSPVL.list)
boxplot(ageinf.list)
boxplot(prev.list)
boxplot(numinc.list)
boxplot(agematch.list); abline(h=0)
boxplot(meanageinf.list)
```
OK, last one was very slightly above, and that too with slight reverse homophily. And with very low incidence.  But it's clear now that things are working, and that the 10 longer runs gives us a much more stable picture of the comparative results.

The time is now ripe to do the deep dive into the treatment. And just overall to see how things are working. Especially the part about how maxcoitalage is about the average age in the relationship, etc. 

Then: consider a run where young people have many shorter relationships and older people have fewer longer relationships.  Right now we're fixing it all with coital frequency, which isn't quire right. And this will perhaps allow for more incidence, and expand the degree to which risk is very different for young than for old.

No wait, that won't work - because that would put the younguns at the far left of the partnership by role graph, but put the old folks further to the right - even if they're on the upper line, they're still below the young folks.  VERY INTERESTING.

OK, deep dive time, then we can return to explorations.

### Dec 13, 2018

Steve went back and read the lit - see table in Evonet/manuscripts/age and SPVL folder.

Discussion  in Evonet meeting:
  - look at other viruses (e.g. HBV, etc)
  - look at whether older people have lower CD4 count, either among HIV-negs or among very newly infecteds (Molly sent paper with hawes and Manhart in the author list)
  - Sarah did quick analysis of MACS data and sawgentle popstive slope in line with Hollingsworth and the like. (but not significant, 0.0047 per year)
  - Josh sent Pantazis paper
  - Josh mentions he and Geoff analyzed MACS data, might be in there.
  
### Mar 04, 2019

Return after a long period of time.

Have been doing literature reviews, brainstorming more components, and attempting to develop a heterosexual model, since most studies seem to be in that population.  

Hetero model is complicated, and requires a fair bit of triangulation from John's scripts and Kathryn's documentation (via Neil's giant spreadsheet of parameter values).

In the meanwhile, the CASCADE study shows an MSM effect (albeit with a much smaller effect size), so will continue trying MSM ones as well.

Moved MSM to a new folder. 

### Mar 12, 2019

Yesterday I created a new branch of Evonet called AgeAndSPVL and added in the basic functionality to have agent heterogeneity in SPVL.  Am going to debug it and test it out now.

But first, need to update text in the .Rmd file to reflect the new locations of different files, both in terms of new organization within this repository, and the fact that the newer larger files (runs 32+) are located outside the repository because they are so large. OK, done.

### Mar 13, 2019

Was able to test and see that the basic functionality appeared to be working. Now I have three full runs going: m43-m45. m43 is exactly the same as m42, and has susceptible_var = 0.So it should be able to test the new functionality while still giving the same results as m42.  Then m44 has suceptibility_var = 0.25, and m45 has 0.5. Results should be available within the hour, conditional on the models actually working! 

Run | 42 | 43 | 44 | 45
---- | ---- | ---- | ----  | ----
min age	| 18 | 18 | 18 | 18
maxage | 55 | 55 | 55 | 55
mean_sqrtage_diff | 1.2 | 1.2 | 1.2  | 1.2
mean_sex_acts_per_day | --- | --- | --- | ---
prob_sex_by_age | T | T | T | T
prob_sex_age_19 | 0.4 | 0.4 | 0.4  | 0.4
max_age_sex | 55 | 55 | 55 | 55
relation_dur | 1000 | 1000 | 1000 | 1000
tx_type | "random" | "random" | "random" | "random"
mean_trtmnt_delay | 0 | 0 | 0 | 0
start_treatment_campaign | 1 | 1 | 1 | 1
proportion_treated | 0.5 | 0.5 | 0.5 | 0.5
testing_model | "interval" | "interval" | "interval"  | "interval"
mean_test_interval_male | 365 | 365 | 365 | 365
prob_tx_dropout | 0.05 | 0.05 | 0.05 | 0.05
susceptibility_var | NA (=0) | *0* | *0.25* | *0.50*

```{r m43-m45}
for (i in 43:45) {

  if(i<10) filler <- "0" else filler <- ""
  load(paste("../AgeAndSPVL_oversize/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  
  popatts[[i]] <- obj$pop
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm
  
  agecoef.list[[i]] <- iSPVL.list[[i]] <- ageinf.list[[i]] <- prev.list[[i]] <- 
    numinc.list[[i]] <- agematch.list[[i]] <- meanageinf.list[[i]] <- vector()
  
  for (j in 1:length(popatts[[i]])) {
    agecoef.list[[i]][j] <- lm(log(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0],10)~
                 popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0])$coef[[2]]
    iSPVL.list[[i]][j] <- mean(log10(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0]),na.rm=TRUE)
    ageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection,na.rm=TRUE)
    prev.list[[i]][j] <- tail(popsumm[[i]][[j]]$prevalence,1)
    numinc.list[[i]][j] <- sum(popatts[[i]][[j]]$Time_Inf>0, na.rm=TRUE)
    agematch.list[[i]][j] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
    meanageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0], na.rm=TRUE)
  }
  rm(obj)
}
```

```{r}
boxplot(agecoef.list); abline(h=0)
boxplot(iSPVL.list)
boxplot(ageinf.list)
boxplot(prev.list)
boxplot(numinc.list)
boxplot(agematch.list); abline(h=0)
boxplot(meanageinf.list)
```

### March 14, 2019

Appears to be working!  At long last!!

I am now going to consider ways to bring up inc/prev, including (1) higher prob_sex_age_; (2) reducing back to durations = 500.


Run | 46 | 47 | 48 | 49 
---- | ---- | ---- | ----  | ----
min age	| 18 | 18 | 18 | 18
maxage | 55 | 55 | 55 | 55
mean_sqrtage_diff | 1.2 | 1.2 | 1.2  | 1.2
mean_sex_acts_per_day | --- | --- | --- | ---
prob_sex_by_age | T | T | T | T
prob_sex_age_19 | 0.4 | 0.4 | 0.8 | 0.8
max_age_sex | 55 | 55 | 55 | 55
relation_dur | 500 | 500 | 500 | 500
tx_type | "random" | "random" | "random" | "random"
mean_trtmnt_delay | 0 | 0 | 0 | 0
start_treatment_campaign | 1 | 1 | 1 | 1
proportion_treated | 0.5 | 0.5 | 0.5 | 0.5
testing_model | "interval" | "interval" | "interval"  | "interval"
mean_test_interval_male | 365 | 365 | 365 | 365
prob_tx_dropout | 0.05 | 0.05 | 0.05 | 0.05
susceptibility_var | 0 | 0.5 | 0 | 0.5

```{r m46-m49}
for (i in 46:49) {

  if(i<10) filler <- "0" else filler <- ""
  load(paste("../AgeAndSPVL_oversize/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  
  popatts[[i]] <- obj$pop
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm
  
  agecoef.list[[i]] <- iSPVL.list[[i]] <- ageinf.list[[i]] <- prev.list[[i]] <- 
    numinc.list[[i]] <- agematch.list[[i]] <- meanageinf.list[[i]] <- vector()
  
  for (j in 1:length(popatts[[i]])) {
    agecoef.list[[i]][j] <- lm(log(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0],10)~
                 popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0])$coef[[2]]
    iSPVL.list[[i]][j] <- mean(log10(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0]),na.rm=TRUE)
    ageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection,na.rm=TRUE)
    prev.list[[i]][j] <- tail(popsumm[[i]][[j]]$prevalence,1)
    numinc.list[[i]][j] <- sum(popatts[[i]][[j]]$Time_Inf>0, na.rm=TRUE)
    agematch.list[[i]][j] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
    meanageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0], na.rm=TRUE)
  }
  rm(obj)
}
```

```{r}
boxplot(agecoef.list); abline(h=0)
boxplot(iSPVL.list)
boxplot(ageinf.list)
boxplot(prev.list)
boxplot(numinc.list)
boxplot(agematch.list); abline(h=0)
boxplot(meanageinf.list)
```

Run | 50 | 51 | 52 | 53
---- | ---- | ---- | ----  | ----
min age	| 18 | 18 | 18 | 18
maxage | 55 | 55 | 55 | 55
mean_sqrtage_diff | 1.2 | 1.2 | 1.2  | **0.7**
mean_sex_acts_per_day | --- | --- | --- | ---
prob_sex_by_age | T | T | T | T
prob_sex_age_19 | 0.6 | 0.6 | **0.8** | 0.6
max_age_sex | 55 | 55 | 55 | 55
relation_dur | 500 | 500 | **1000** | 500
tx_type | "random" | "random" | "random" | "random"
mean_trtmnt_delay | 0 | 0 | 0 | 0
start_treatment_campaign | 1 | 1 | 1 | 1
proportion_treated | 0.5 | 0.5 | 0.5 | 0.5
testing_model | "interval" | "interval" | "interval"  | "interval"
mean_test_interval_male | 365 | 365 | 365 | 365
prob_tx_dropout | 0.05 | 0.05 | 0.05 | 0.05
susceptibility_var | 0.5 | **1** | **1** | **1**

```{r m50-m53}
for (i in 50:53) {

  if(i<10) filler <- "0" else filler <- ""
  load(paste("../AgeAndSPVL_oversize/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  
  popatts[[i]] <- obj$pop
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm
  
  agecoef.list[[i]] <- iSPVL.list[[i]] <- ageinf.list[[i]] <- prev.list[[i]] <- 
    numinc.list[[i]] <- agematch.list[[i]] <- meanageinf.list[[i]] <- vector()
  
  for (j in 1:length(popatts[[i]])) {
    agecoef.list[[i]][j] <- lm(log(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0],10)~
                 popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0])$coef[[2]]
    iSPVL.list[[i]][j] <- mean(log10(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0]),na.rm=TRUE)
    ageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection,na.rm=TRUE)
    prev.list[[i]][j] <- tail(popsumm[[i]][[j]]$prevalence,1)
    numinc.list[[i]][j] <- sum(popatts[[i]][[j]]$Time_Inf>0, na.rm=TRUE)
    agematch.list[[i]][j] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
    meanageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0], na.rm=TRUE)
  }
  rm(obj)
}
```

```{r}
boxplot(agecoef.list); abline(h=0)
boxplot(iSPVL.list)
boxplot(ageinf.list)
boxplot(prev.list)
boxplot(numinc.list)
boxplot(agematch.list); abline(h=0)
boxplot(meanageinf.list)
```

Summarizing the less few sets of runs, it would seem that:

- susceptible_var is indeed important
- age mixing is indeed important
- prob of sex is not that important for age coefficient (but hugely impt for prevalence)
- relational duration is not that important for age coefficient (but hugely impt for prevalence)

Next steps:
- Add in susceptibility attribute in remaining places (bookkeeping stuff)
- Review whole annals so far and make list of all loose ends / ideas / etc.
- Check on the treatment levels and cycling to confirm the code is behaving as understood
- Consider upper limit of susceptible_var reasonable
- check on when people age out


- Begin to set up multi-dimensional experiments on 10k MSM pops?
- Pick back up on het pops?

### April 19, 2019

Back to it. First off, i'm trying to bring in the changes to evonet that James made to allow it to work without fast_edgelist. I do so with 

git checkout AgeAndSPVL
git merge origin/master

my understanding is that this brings the changes from origin/master into AgeAndSPVL, and the commit history would seem to confirm that.

Oddly, though, both branches on my local machine are getting caught on the line export(treatment_dropout_test) in NAMESPACE.  When I comment it out, things work.  I need to determine what exactly is going on here; from the name it sounds like something I would have added to the AgeAndSPVL branch, but why then would it show up in master branch as well?  And when I changed it in one, how come it is showing up in the git window as a pending change regardless of which branch I toggle to?

In the meanwhile, I am testing whether the code works as expected (modulo the above issue, which I commented out for now).  I am doing so by creating a new run (m54) which is just m27 but with the births and deaths modules replaced with evo_arrivals and evo_departures, as James told us to do. I can them look to see if the results are of the same flavor.  If this works (I can go through and change all scripts to do this; or perhaps just change them moving forward actually.

### April 24, 2019

The treatment_dropout_test issue originated with James; he fixed it, and I just updated Evonet (master and AgeAndSPVL branches).  SO now all should be back on track, and I can start confirming that treatment dropout occurs as I expect.  (The overlap between my interest in testing treatment drop and James's is pure coincidence).

### April 30, 2019

Indeed, run m54 does look just like run 27, so that is all good.

Then spent some time a few days ago confirming that treatment dropout does occur as expected, and the answer is that it does.

Now going to try to implement treatment initiation using a fixed percent rule, following the code that John laid out in his email dated 4/10.

Have added the percent_of_current_diagnosed option for tx_limit.  It is now running, as m55 (identical to m31 - the one I analyzed and discussed in the email on 4/10 - but with percent_of_current_diagnosed). We shall see!

### May 1, 2019

run 55 finished, and it appears to behave exactly as expected, with: 

```
aaa <- ageSPVL_m55$popsumm[[1]]
plot(aaa$no_treated / aaa$total_infections_alive)
```

and when I look at the vl plots, and when i look at the popsumm figures.  

Now I have to decide about an age distribution to start with.  i'l lbe varying treatment levels, so there's no way to make it work with all of them. I just need to decide what is a reasonable place overall.  I'm thinking just using ASMR.

OK, found the USA-Male numbers, and tried to get that to work, but had a fair bit of trouble.  Outlined in an email I sent to theteam on 5.1/  John wrote back with some code that worked for him, and an explanation.  iI implemented that in run 56, and it worked as advertised.  so i am going with that code for now.

Setting up run 57 now to test the new functionality on a set of parameters with reasonable likelihood of a positive correlation and sizeable epidemic.  

Also setting up run 58 to match 57 but have a max date of 75 in it. Sarah had suggested in a meeting that maybe that would help - even though max age sex is 55, that only refers to the average age.

Run | 57 | 58
---- | ----
min age	| 18 | 18
maxage | 55 | 75
mean_sqrtage_diff | 0.7 | 0.7
prob_sex_age_19 | 0.8 | 0.8
max_age_sex | 55 | 55 
relation_dur | 1000 | 1000
tx_type | "random" | "random"
mean_trtmnt_delay | 0 | 0 
start_treatment_campaign | 3 | 3
proportion_treated | 0.5 | 0.5 
testing_model | "interval" | "interval"
mean_test_interval_male | 365 | 36
prob_tx_dropout | 0.05 | 0.05
susceptibility_var | 1.0 | 1.0
tx_limit | "percent_of_currently_diagnosed" | "percent_of_currently_diagnosed"
initial_agedata_male | "stable_age_no_hiv_dist" | "stable_age_no_hiv_dist"
birth_model | "exponential_growth" | "exponential_growth"
baseline_input_exp_growth | 0.007*(evoparams$initial_pop/100) | 0.007*(evoparams$initial_pop/100) 


```{r m57-m58}
for (i in 57:58) {

  if(i<10) filler <- "0" else filler <- ""
  load(paste("../AgeAndSPVL_oversize/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  
  popatts[[i]] <- obj$pop
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm
  
  agecoef.list[[i]] <- iSPVL.list[[i]] <- ageinf.list[[i]] <- prev.list[[i]] <- 
    numinc.list[[i]] <- agematch.list[[i]] <- meanageinf.list[[i]] <- vector()
  
  for (j in 1:length(popatts[[i]])) {
    agecoef.list[[i]][j] <- lm(log(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0],10)~
                 popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0])$coef[[2]]
    iSPVL.list[[i]][j] <- mean(log10(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0]),na.rm=TRUE)
    ageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection,na.rm=TRUE)
    prev.list[[i]][j] <- tail(popsumm[[i]][[j]]$prevalence,1)
    numinc.list[[i]][j] <- sum(popatts[[i]][[j]]$Time_Inf>0, na.rm=TRUE)
    agematch.list[[i]][j] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
    meanageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0], na.rm=TRUE)
  }
  rm(obj)
}
```

```{r}
boxplot(agecoef.list); abline(h=0)
boxplot(iSPVL.list)
boxplot(ageinf.list)
boxplot(prev.list)
boxplot(numinc.list)
boxplot(agematch.list); abline(h=0)
boxplot(meanageinf.list)
```
### May 3, 2019

Greetings from Maine. Ok, runs 57 and 58 worked, but produced bland results.  Of course I did my usual: I was so excited that I was on the track of something interesting that I changed multiple things relative to any previous run, so now it is impossible to isolate effects.  

I am going now going to do four runs:

- 59: matches 52 exactly, but will be a test to make sure that the changes to the code didn't create any unexpected differences
- 60: 59, plus the new tx_limit
- 61: 60, plus the new init_age and associated arguments (four in total!!!!)
- 62: 61, plus moving max age to 75.

Run | 59 | 60 | 61 | 62
---- | ---- | ---- | ----  | ----
min age	| 18 | 18 | 18 | 18
maxage | 55 | 55 | 55 | 75
mean_sqrtage_diff | 1.2 | 1.2 | 1.2  | 1.2
mean_sex_acts_per_day | --- | --- | --- | ---
prob_sex_by_age | T | T | T | T
prob_sex_age_19 | 0.8 | 0.8 | 0.8 | 0.8
max_age_sex | 55 | 55 | 55 | 55
relation_dur | 1000 | 1000 | 1000 | 1000
tx_type | "random" | "random" | "random" | "random"
mean_trtmnt_delay | 0 | 0 | 0 | 0
start_treatment_campaign | 3 | 3 | 3 | 3
proportion_treated | 0.5 | 0.5 | 0.5 | 0.5
testing_model | "interval" | "interval" | "interval"  | "interval"
mean_test_interval_male | 365 | 365 | 365 | 365
prob_tx_dropout | 0.05 | 0.05 | 0.05 | 0.05
susceptibility_var | 1.0 | 1.0 | 1.0 | 1.0
tx_limit | "absolute_num" | "percent_of_currently_diagnosed" | "percent_of_currently_diagnosed" | "percent_of_currently_diagnosed"
initial_agedata_male | "linear_decrease" | "linear_decrease" | "stable_age_no_hiv_dist" | "stable_age_no_hiv_dist"
birth_model | default | default | "exponential_growth" | "exponential_growth"
baseline_input_exp_growth | default | default | 0.007*(evoparams$initial_pop/100) | ditto  


```{r m59-m62}
for (i in 59:62) {

  if(i<10) filler <- "0" else filler <- ""
  load(paste("../AgeAndSPVL_oversize/ageSPVL_m",filler,i,".rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,sep=""))
  
  popatts[[i]] <- obj$pop
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,sep=""))$popsumm
  
  agecoef.list[[i]] <- iSPVL.list[[i]] <- ageinf.list[[i]] <- prev.list[[i]] <- 
    numinc.list[[i]] <- agematch.list[[i]] <- meanageinf.list[[i]] <- vector()
  
  for (j in 1:length(popatts[[i]])) {
    agecoef.list[[i]][j] <- lm(log(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0],10)~
                 popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0])$coef[[2]]
    iSPVL.list[[i]][j] <- mean(log10(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0]),na.rm=TRUE)
    ageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection,na.rm=TRUE)
    prev.list[[i]][j] <- tail(popsumm[[i]][[j]]$prevalence,1)
    numinc.list[[i]][j] <- sum(popatts[[i]][[j]]$Time_Inf>0, na.rm=TRUE)
    agematch.list[[i]][j] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
    meanageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0], na.rm=TRUE)
  }
  rm(obj)
}
```

```{r}
boxplot(agecoef.list); abline(h=0)
boxplot(iSPVL.list)
boxplot(ageinf.list)
boxplot(prev.list)
boxplot(numinc.list)
boxplot(agematch.list); abline(h=0)
boxplot(meanageinf.list)
```

### May 14, 2019

Runs 59-62 showed nice results.  Especially the fact that, by including up to age 75, I am able to bring the agediff coef down below 0, while keeping the actual coital acts by age pretty much the same, and thus the age/SPVL correlation results.  Nice!

61 was the highest, followed closely by 62. So that means John's code to get the age distribution right seems to have helped a lot. And even though 62 is lower than 61, 62 is the run with the negative absdiff term.  Indeed I believe that is the fist ever negative absdiff term to have a postive age/SPVL correlation. Which is a great step forward.

I am now going to send Sarah a clean starting script and sets of parameters to explore.

It will still be worth thinking about additional ways to vary risk by age, probably with Johns' heterosexual code.  E.g. shorter relationships and higher concurrency for younger folks.

OK, I have created a list of runs for Sarah which I will email to her now.  It is "_supporting_files/AgeSPVL runs master list.xlsx".   I've also created the first of the script files (m101).

### May 23, 2019

Sarah has done the runs.

Instructions to get the results is:

0. One-time: make directory in libra called evonet/AgeAndSPVL/runs
1. Open terminal window in RStudio.
2. Type ssh -X goodreau@mox.hyak.uw.edu
3. Enter password, and follow Duo instructions.
4. One time instructions: set up .ssh/config file, as per instructions at https://github.com/statnet/computing/wiki/Getting-Started-with-Hyak#ssh-config-file and in Sarah's email from 5/22/2019
5. Type scp /gscratch/csde/sestans/age/* libra:~/evonet/AgeAndSPVL/runs
6. All runs should be on libra

Now beginning analysis on libra, by copying over the above analysis code there.  

### May 28, 2019

Did the initial analyses of runs m101-m149, with a few runs still missing.  But overall the patterns were very interesting.  Most things behaved as expected.  Percent on treatment had a big effect, although rate going off did not.  The weirdest part was how little impact there was in the exploration of susceptibility_var.  indeed, the run with susceptibility_var of 0 was still a positive effect, not something I've ever seen efore.  This requires deeper exploration ASAP. I shared these in our Evnet call on Friday, with some good feedback (see todo list below).

Sarah has set more replicates going of all of these scenarios, so that we should soon go from 24 to 96.  I am going to compile those that have finished now, to help me see the patterns more clearly.  Then there will be some things to try.

### May 29, 2019

On my way to reunion; compiled all reulsts for m101-149. Delved deeper into the weird susceptibility_var results.  Turns out that the way I implemented means that even the values I was looking at (e.g. var = 1 and above) had enormous amouts of variance. It's in part because I called it var but then impemented as sd, but then more importantly because in the Hughes formula the then get exponentiated. I need to dive back into that.  But first I am rushing to create a new set of runs for Sarah to produce while I'm at the reunion - m150-191.  These are mostly all the same scenarios as m101-m149 but with no susvar = 0 instead of 1.  Then there is a few that re-explore susvar values but at a more fine-grained level between 0 and 1.

### Jun 4, 2019

Compiled all of the runs from the reunion days (m150-m191). The patterns are remarkably similar to the previous set that had susvar=1. This means there are two inter-related questions about the susvar hanging out there:
- why does it not seem to matter when earlier on during my epxlorations I thought it did?
- why do the runs that have susvar=0 seem to almost all have positive coefficients now, when before I thought they had negative?

These questions will need to be dived into again before publication - but for now (upcoming Sunbelt paper, upcoming Epidemics abstract deadline), I can put them aside.

more important is to consider a few more network cinsiderations.  First and foremost is shorter relations for youngins.

Creating run 192, which matches run 150, but with split durations. Setting up syntax to try locally, then if it looks promisin, will send to Sarah.


### July 12, 2019

Back from crazy month.  Step 1 is to ascertain why the run 192 code worked for me and not for Sarah.  James had identified that it worked for him with fastedgelist+F but not with T.

Steve and Sarah compared sessionInfos, Steve discovered quite a few differences, and then he stumbled on a file called AgeAndSPVL_package_installs.R that had the lines

remotes::install_github("https://github.com/chad-klumb/network", ref="tergmLiteterms")
remotes::install_github("https://github.com/statnet/tergmLite", ref="tergmLiteterms")
remotes::install_github("https://github.com/statnet/EpiModel", ref="tergmLiteterms")

thus triggering his memory that these are needed to make the new functionality work.  he installed some and then it stopped working.  Now trying to figure out the right combination.....



### Aug 29, 2019

CAMP is under control again, so returning here.

Right now I am stepping back and taking stock of all of the runs and what they found, before moving forward with the final push.  A quick review last week suggested that I over-estiamted the impact that sus_var had - it really might not have had much of all.  Now I am looking into see howif I can create some of the highest parameter runs from my local machine onto Hyak, and comparing them with the most similar runs already done at Hyak.

Specically I am focusing on run 43 on the local side, and run 151 on the Hyak side.  From my initial review of parameters, they appear to differ only in the probability of sex at age 19 (0.4 in run 43, 0.8 in run 151).  Needless to say, run 43 has a much smaller prevalence and incidence than run 151 does.  But it also has a higher regression  coefficient, at least in terms of median (0.0053, although mean is lower at 0.0028, with 10 runs). This is compared to run 151, with 0.0018 for both mean and median).  So I am going to try now to run run 43 as is on Hyak, then run a new model beefing up size and time, then a third model that just changes run 43 to have 0.8 prob sex at 19.

### Aug 30, 2019

All meaningful parameter differences between m43 and m151:

parameter | m43 | m151
-------- | --- | -----
hyak args | no | yes 
initial_pop | 1000 | 10000
prob_sex_age_19 | 0.4 | 0.8
initial_agedata_male | "linear_decrease" | "stable_age_no_hiv_dist"
initial_infected | 200 (20%) | 1000 (10%)
tx_limit | "absolute_num" | "percent_of_currently_diagnosed"
  
I am thus creating new runs now, changing one at a time:
 
parameter | m43 | m193 | m194 | m195 | m196 | m197 | m151 
-------- | --- | --- | --- | ---  | --- | --- | ---
hyak args | no | *yes* | yes | yes | yes | yes | yes
initial_agedata_male | "linear_decrease" | "linear_decrease" | *"stable_age_no_hiv_dist"* | "stable_age_no_hiv_dist" | "stable_age_no_hiv_dist" | "stable_age_no_hiv_dist" | "stable_age_no_hiv_dist"
tx_limit | "absolute_num" | "absolute_num" | absolute_num" | *"percent_of_currently_diagnosed"* | "percent_of_currently_diagnosed" | "percent_of_currently_diagnosed" | "percent_of_currently_diagnosed" 
initial_infected | 200 (20%) | 200 (20%) | 200 (20%) | 200 (20%) | *100 (10%)* | *1000 (10%)* | 1000 (10%)
initial_pop | 1000 | 1000 | 1000 | 1000 | 1000 | *10000* | 10000
prob_sex_age_19 | 0.4 | 0.4 | 0.4 | 0.4 | 0.4 | 0.4 | *0.8*

Note that initial_infected has to be changed twice: once to represent the change in intial prevalence ,and once to represent the change in initial population size.


### Sep 9, 2019

Results of these runs are promising.  High coefficients, although wide CIs since they have the small pop sizes.  I am going to rerun them now with n=10k and prob_sex_age_19 = 0.8 to see.  Then the time has come to look into run 192 again.

parameter | m198 | m199 | m200 | m201 
-------- | --- | --- | ---  | --- 
hyak args | yes | yes | yes | yes
initial_agedata_male | "linear_decrease" | *"stable_age_no_hiv_dist"* | "stable_age_no_hiv_dist" | "stable_age_no_hiv_dist"
tx_limit | "absolute_num" | absolute_num" | *"percent_of_currently_diagnosed"* | "percent_of_currently_diagnosed"
initial_infected | 20% | 20% | 20% | *10%*
initial_pop | 10000 | 10000 | 10000 | 10000
prob_sex_age_19 | 0.8 | 0.8 | 0.8 | 0.8

Remember that these all have maxage 55, with age diff 1.2.  This is all part of comparing runs between libra and hyak, especially those with strong effects.  But will need to increase maxage much higher, and add in the rapidly declining coital frequency, once these are done.

### Sep 19

Last week I hit a wall where suddenly my code was not working on Hyak.  Along the way I learned that the `tergmLiteterms` branch of `EpiModel` had been merged into `master`, necessitating one change in my set-up instructions:

remotes::install_github("https://github.com/chad-klumb/network", ref="tergmLiteterms")
remotes::install_github("https://github.com/statnet/tergmLite", ref="tergmLiteterms")
remotes::install_github("https://github.com/statnet/EpiModel", ref="master")

However, updating that led to an error about "sqrt_age" not being a nodal attribute.  Sarah informed me she had gotten the same error once, and reverting to old versions of packages had solved it. I could try that as a stopgap (although getting all the right versions that can handle the "~age<25" argument is tricky; plus it doesn't solve the problem.  So now I am going to try to create stand-aone-is scrips that show the problem on my laptop, where I can get it to work with one old version (without age<25) but not with the newest versions.

Sarah emailed versions of the packages (already unpacked, not as tarballs) that she says works, and that match my memory of what should work.  They are currently in C:\active\EvoNet\manuscripts\age and SPVL\temp r pkgs\EpiModel.  Based o ntheir DESCRIPTION files, they are:

Package: EpiModel
Version: 1.7.3
Date: 2019-04-24

Package: network
Version: 1.16-377
Date: 2019-04-01

Package: tergmLite
Version: 1.2.0
Date: 2018-11-02

I have created a version of m193a to run locally called m193_localdebug.  This file worked just a few weeks ago on Hyak. I wil lsee if I can get it to run with either the versions of the packages I have here now or soemthing similar, and then to fail with the most recent versions.  Then I will see if I can pair the whole thing down to relatively little code.

### Sep 25, 2019

Lots of sleuthing by Chad has found the bug - something related to Roxygen and exporting generics and Sam not having Chad's version of network on his computer when he ran Rxoygen for tergmLite.Really weird stuff. But it's fixed now.

So the time has come to send runs 198-201 to Hyak again.

All look good.  High coefs, small intervals.  So now we can finally finally move forward again.

From here on out, will use: 

initial_agedata_male = *"stable_age_no_hiv_dist"*
tx_limit = *"percent_of_currently_diagnosed"*

And will return to higher maxage in population with still lower maxage of sex.

### Oct 14, 2019

Have spent time on and off over the last few weeks looking into 192 more deeply.  Code does indeed seem to be working.  Parameterizing the models are tricky though; one needs to have nodemix(age>25) in the formation to go along with it in the dissolution.  And ocming up with reasonable target stats when the model also has absdiff(age) is tricky.  For now pausing on that and going into non-linear declines in coital frequency by age.

Branch for trying this added to EvoNetHIV (coitalfreq_additions).  Code written, using a new options:

prob_sex_func: TRUE
coital_coef_base: logodds of sex per time step for two 18-year olds in a new relationship
coital_coef_age  = decrease in logodds per year of age for either member
coital_coef_dur  = decrease in logodds per year of relational duration (not yet implemented)

Tried out an exmaple on the sim cluster (ageSPVL_m202_local.R) with values of coital_coef_base = 1.4 and coital_coef_age = -0.6.  seems to work well, with a coef of 0.002.  Will try again now with ageSPVL_m203_local.R and coital_coef_age = -0.8.  If results are sensible I will upgrade to full runs on Hyak.


### Oct 22, 2019

Sim 203 gave unexciting results.  I then rewrote the prob_sex_func code to take a different form: the lower-half of a declining logistic, with both a multiplier and an additive applied after the converstion to a probability.  THis allows for arbitrary stating probability and also arbitrary asymptote.  I first tried it with a small asymptote (5%? m204) and, with disappointing results, did it again with an asymptote of 0 (again as m204).  Results were similar to 203 - slightly under 0.001.  This contrasts substantially with run m168 (max age sex = 45, age coef = 0.003), even though they should be rights be pretty similar.  So I will now look a bit closer at both.

In the meanwhile, I was reminder that in an Evonet meeting that (1) setting the flag to allow the acute peak to covary with SPVL might accentuate my results, and (2) increasing the acute stage multiplier might up my incidence and prevalence.

Only differences in m204 vs m168 were: (1) use of hyak for latter; (2) the different age coital code of interest; (3) m204 didn't set min_age, although the valuewe use--18--is the default, so it shouldn't matter, and (4) initial prevalence was 10% in m168 and 20% in m204.  Just to be sure that there are no meaningful differences, I'm adapting m204 into m205, and setting back the coital code to use the sex_by_age argument.

### Oct 23, 2019

OK, m205 showed the original strong effect.  So either there's a bug in my code or logic for prob_sex_func, or that small difference in coital freq by age actually matters a lot.  I'm guessing the former.  Will dive in now.

### Oct 31, 2019

Happy Halloween.  Just back from State College, and can update this journal with progress over that trip.  

So, it turns out that the exact function of sex by age ends up mattering *a lot*.  I have settled on a version that declines logistically to an asymptote, based on the *harmonic mean* of the ages.  This makes a big difference,  It implies that a relationship between an old and young person is a bit more driven by the young person -- perhaps because it is reasonably likely to be short or novel, or otherwise odd.  Got the biggest efect we've yet seen, as run 206.

Noting, however, that the prevalences here are very low, and also that acute infection differences might help us out, I next decided to explore (1) making the acute peak a function of SPVL, and bumping up the acute peak multiplier.  I did these in runs 207 (1); 208 (1, 2); 209 (2).  All of them saw reductions away form the larger effect in run 206 (although 209 was the least reduced, suggesting that (1) is the main cause of the reductions. Runs 208 and 209 did see increases in prevalence, but only slightly -- I think the relaitonal durations really need to come down to bring that up.

Next I have two directions I can go in: (1) adding in the code to make coital freq a function of relational age; and (2) re-exploring models where rel dur is shorter for younger folks.  I am going to do the latter first because I know it doesn't  require any mucking around with the inner workings of evonet -- I'd already gotten it to work before, although with model parameters that were unrealistic. 

The first part of this will be in revisiting exactly how the current pattern of age mixing in the model (absdiff sqrt.age = 1.2) behaves given the initial age distribution we've settled on.  It yields a negative coefficient, but not an extremely strong one -- but plots of age mixing taken during the debug phase -- at least for the discordant edgelist at the start of the model, which should be close to a random selection of relationships -- show a pretty scattered relational distribution.  A little more sense of this will be needed before plunging into changing both the formation and dissolution rates by age. I should also take a look at mean degree in the process.

Important note:  since all of this work is being done in an exploratory fashion here on the sim cluster instead of on Hyak -- I am writing the files in a speciak folder, with sepcial names (originally containing "debug" and now containing "local").  The output is too large to save into git, so I am copying it out to Age_and_SPVL_oversize.

OK, dove in to see what sqrt age 1.2 really means.  Did so with run 206 and the debugger.  First off, it runs as expected - absdiff sqrt age in the initial network of the sim is 1.19.  It turns out that this corresponds to an absolute mean age difference of 15.1 !!  So yeah. that's a big age difference.  Hence the wide scatter.  We will stick with it for some time to see how high we can go, but then we will need to dial it back as well during the final runs.

So now comes the time to settle on terms for the model with declining durations by age.

Q1: can I use a nodecov in the dissolution and formation and make my life easier?  I think that I may be able to; I just can't use the built-in dissolution_coef function.  But that should be fine, no? Exploring now.

### Nov 6, 2019

OK, so run 210 (with nodemix on formation and dissolution for age>25) yields a coef close to 0. I have the model right so that the nodexis coefs prior to the Carnegie correction are basically 0, meaning the density is as expected if thereall groups have the same mean degree.  But then the mean degree does drift off after time passes - md for under 30 is right on 0.7, but for 30-50 it rises to 0.8, and for 50+ it rises to 0.75.  Interesting.

Meanwhile, I did a review of the earlier models to see where if anywhere I could get a bit more leverage. Switching from duration 1000 to 500 might get more prevalence without reducing the effect size much or at all.  And raising the dropoff rate to 0.2. And raising mean degree.

I am going to try runs for each of those now, as such:

211: youth have duration 500, older 1000, mixed = harmonic mean = 667.
212: 206 but with dropout = 0.2
213: 212 but with mean deg = 0.8
214: 213 but with relation dur = 500

Then I will look to see about adding in nodecov(age) instead of durations on age.


### Nov 11, 2019

Ok, none of those were exciting.  Indeed, 212 pretty much caused the effect to go away, which is very weird.

215: 206 but with meandeg = 0.8
216: 206 but with meandeg = 0.9
217: 206 but with nodecov(age)
218: 206 but with concurrency at 5%  - never thought to try this yet !!
219: 206 but with concurrency higher for younger people
220: 206 but only 20 years long


```{r m202-218}
for (i in 202:218) {

  if(i<10) filler <- "0" else filler <- ""
  load(paste("../AgeAndSPVL_oversize/ageSPVL_m",filler,i,"_testa.rda",sep=""))
  obj <- get(paste("ageSPVL_m",filler,i,'_testa',sep=""))
  
  popatts[[i]] <- obj$pop
  popsumm[[i]] <- get(paste("ageSPVL_m",filler,i,'_testa',sep=""))$popsumm
  
  agecoef.list[[i]] <- iSPVL.list[[i]] <- ageinf.list[[i]] <- prev.list[[i]] <- 
    numinc.list[[i]] <- agematch.list[[i]] <- meanageinf.list[[i]] <- vector()
  
  for (j in 1:length(popatts[[i]])) {
    agecoef.list[[i]][j] <- lm(log(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0],10)~
                 popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0])$coef[[2]]
    iSPVL.list[[i]][j] <- mean(log10(popatts[[i]][[j]]$SetPoint[popatts[[i]][[j]]$Time_Inf>0]),na.rm=TRUE)
    ageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection,na.rm=TRUE)
    prev.list[[i]][j] <- tail(popsumm[[i]][[j]]$prevalence,1)
    numinc.list[[i]][j] <- sum(popatts[[i]][[j]]$Time_Inf>0, na.rm=TRUE)
    agematch.list[[i]][j] <- obj$nwparam[[1]]$coef.form['absdiff.sqrt_age']
    meanageinf.list[[i]][j] <- mean(popatts[[i]][[j]]$age_infection[popatts[[i]][[j]]$Time_Inf>0], na.rm=TRUE)
  }
  rm(obj)
}
```


```{r}

xlim <- c(201,220)
boxplot(agecoef.list, xlim=xlim); abline(h=0)
boxplot(iSPVL.list, xlim=xlim)
boxplot(ageinf.list, xlim=xlim)
boxplot(prev.list, xlim=xlim)
boxplot(numinc.list, xlim=xlim)
boxplot(agematch.list, xlim=xlim); abline(h=0)
boxplot(meanageinf.list, xlim=xlim)
```

### May 7, 2020

Picking up once again after a long long time (thanks for nothing, COVID).

We are moving into the "fish or cut bait" phase of this paper.

Looks like I did a whole ton of runs in the 300 range, and did not document them here.I think this was for the Epidemics talk, but will figure out now.


### To do (always keep at the bottom)

To Do:
- consider decline with rel length 
- add shorter relational durations for younger folks.
- add more concurrency for younger folks?
- #heteros!
- #revisit susvar
- consider acute flag and multiplier.