R/summary_popsumm.R

In EvoNetHIV/Herbeck-et-al-Vaccine-201x: Evonet: Integrated Bio-Social Models for HIV Epidemiology

#' @export
summary_popsumm<-function(dat,at){

  
  #if summary stats calculations doesn't occur at this timesetp,
  #based on popsumm_frequency value, then leave function
  if( (at%%dat$param$popsumm_frequency!=0) & (at!=1)){return(dat)}
  
  #time_index is a time vector based on current value
  #of "at" and parameter value "popsumm_frequency"
  if(at==1)
    time_index <- 1
  else if(at==dat$param$popsumm_frequency)
    time_index <- (at-dat$param$popsumm_frequency+2):at
  else
    time_index<- (at-dat$param$popsumm_frequency+1):at
  
  #"popsumm_index" is an index for the "popsumm" vectors
  #based on value of "at" and paramter "popsumm_frequency"
  if(at==1)
  popsumm_index <- 1
  else
    if(dat$param$popsumm_frequency==1)
      popsumm_index <- at
    else
      popsumm_index <- (at/dat$param$popsumm_frequency)+1
  
  
  #logical vectors and indices helpful to calculate summary stats  
    
  inf_index     <-  dat$pop$Status == 1
  total_inf     <- length(which(inf_index))
  sus_index     <-  dat$pop$Status == 0
  alive_index   <-  inf_index | sus_index
  total_alive <- length(which(alive_index))
  treated_index <-  dat$pop$treated == 1 & inf_index
  not_treated_index <-  dat$pop$treated == 0 & inf_index
  treated_undetectable <- treated_index & dat$pop$V<dat$param$vl_undetectable
  treated_agents <- which(treated_index)
  not_treated_agents <- which(not_treated_index)
  
  under30_index <- dat$pop$age < 30
  agents30to50_index <- dat$pop$age >= 30 & dat$pop$age < 50 # Testing need to change 31 back to 50
  over50_index <- dat$pop$age >= 50
  agents_under30 <- which(under30_index)
  agents_30to50 <- which(agents30to50_index)
  agents_over50 <- which(over50_index)
  
  new_infections <- is.element(dat$pop$Time_Inf, time_index)
  new_infections_count <- length(which(is.element(dat$pop$Time_Inf, time_index)))
  new_infections_virus_vacc_sens_count <- length(which(is.element(dat$pop$Time_Inf, time_index)&
                                                         dat$pop$virus_sens_vacc==1))
  new_infections_virus_vacc_notsens_count <- length(which(is.element(dat$pop$Time_Inf, time_index)&
                                                         dat$pop$virus_sens_vacc==0))
  
  new_infections_virus_drug_sens_count <- length(which(is.element(dat$pop$Time_Inf, time_index)&
                                                         dat$pop$virus_sens_drug==1))
  new_infections_virus_drug_part_res_count <- length(which(is.element(dat$pop$Time_Inf, time_index)&
                                                         dat$pop$virus_part_res_drug==1))
  new_infections_virus_drug_3_plus_res_count <- length(which(is.element(dat$pop$Time_Inf, time_index) &
                                                         dat$pop$virus_3_plus_drug_muts==1))
  new_infections_virus_1_drug_muts <- length(which(is.element(dat$pop$Time_Inf, time_index) &
                                                             dat$pop$virus_3_plus_drug_muts==1))
  
  
  donor_time_inf  <- ifelse(new_infections_count>0,
                            dat$pop$Donors_Total_Time_Inf_At_Trans[new_infections],
                            NA)
  donor_acute_count <- ifelse(!is.na(donor_time_inf),
                              length(which(donor_time_inf<=dat$param$t_acute)),
                              NA)
  new_births <- is.element(dat$pop$arrival_time, time_index)
  cd4_aids <- dat$pop$CD4 == 4
  new_diagnoses <- dat$pop$diag_status == 1 &  is.element(dat$pop$diag_time,time_index)
  acute_phase_vec <- (at-dat$pop$Time_Inf)<dat$param$t_acute
  acute_phase <- !is.na(acute_phase_vec) & acute_phase_vec==T
  percent_virus_sensitive <- round(100*(length(which(dat$pop$virus_sens_vacc==1 & inf_index))/length(which(inf_index))))
  percentVaccinated <- round(100*(length(which(dat$pop$vaccinated == 1 & alive_index))/total_alive))
  
  #deaths
  just_died <- is.element(dat$pop$Time_Death,time_index)
  died_aids <- dat$pop$Status == -2 & just_died
  died_aids_mean_age <- mean(dat$pop$age[died_aids])
  died_non_aids <- dat$pop$Status == -1 & just_died
  died_non_aids_inf <- died_non_aids & !is.na(dat$pop$V)
  died_non_aids_sus <- died_non_aids & is.na(dat$pop$V)
  aged_out <- (dat$pop$age>=dat$param$max_age) & just_died
  
#browser()
  #network statistics
  # some of these can't be computed if we are in edgelist mode
  # so need to create a network from the edgelist
  if(!is.null(dat[['nw']])){
    nw <- dat[['nw']]
  } else {
    nw_summary    <-  NULL 
    number_edges <- nrow(dat$el) 
 	  network_size <- attr(dat$el,'n') 
    total_nodes   <-  NULL 
    netattrs<-attributes(dat$el)
    nw <- as.network.matrix(dat$el, matrix.type='edgelist',
                            # TODO: ASSUMING THESE HAVE BEEN HARDCODED UPSTREAM
                            directed = FALSE,
                            bipartite = FALSE,
                            loops = FALSE
                            )
  }
  nw_summary    <-  summary(nw~degree(0:1) + concurrent, at = at)
  number_edges <- network.edgecount(nw)
  network_size <- network.size(nw)
  total_nodes   <-  sum(nw_summary[1]+nw_summary[2]+nw_summary[3]) # This depends on nw_summary which I blanked out above
    
 
  #viral load values
  log10_vl_values  <-  log10(dat$pop$V[which(inf_index)]+dat$param$AbsoluteCut)
  spvl_untreated_values <- (
        dat$pop$LogSetPoint[which(inf_index & not_treated_index)])
  # todo: may be a faster way to calculate degree
  
  edges_by_agent <- unname(summary(nw ~ sociality(base = 0),at=at)) #use dat$attr$id for index on dat$pop
  edges_untreated <- edges_by_agent[dat$attr$id %in% not_treated_agents ]
  edges_treated <- edges_by_agent[dat$attr$id %in%  treated_agents]
  edges_under30 <- edges_by_agent[dat$attr$id %in% agents_under30]
  edges_30to50 <- edges_by_agent[dat$attr$id %in% agents_30to50]
  edges_over50 <- edges_by_agent[dat$attr$id %in% agents_over50]
  
   #aim3 mutations
   inf_undetect_ix <- (dat$pop$Status==1 & dat$pop$V> dat$param$vl_undetectable)
   no_inf_undect <- length(which(inf_undetect_ix))
  mutations0 <- length(which(inf_undetect_ix & dat$pop$aim3_no_muts==0))
  mutations1 <- length(which(inf_undetect_ix & dat$pop$aim3_no_muts>=1))
  mutations2 <- length(which(inf_undetect_ix & dat$pop$aim3_no_muts>=2))
  mutations3 <- length(which(inf_undetect_ix & dat$pop$aim3_no_muts>=3))
  mutations4 <- length(which(inf_undetect_ix & dat$pop$aim3_no_muts>=4))
  mutations5 <- length(which(inf_undetect_ix & dat$pop$aim3_no_muts>=5))
  
  mutations1exact <- length(which(inf_undetect_ix & dat$pop$aim3_no_muts==1))
  mutations2exact <- length(which(inf_undetect_ix & dat$pop$aim3_no_muts==2))
  mutations3exact <- length(which(inf_undetect_ix & dat$pop$aim3_no_muts==3))
  mutations4exact <- length(which(inf_undetect_ix & dat$pop$aim3_no_muts==4))
  
  mutations3plus_long <- length(which(inf_index & dat$pop$aim3_muations_long>=3)) 
  mutations4plus_long <- length(which(inf_index & dat$pop$aim3_muations_long>=4)) 
  mutations5_long <- length(which(inf_index & dat$pop$aim3_muations_long==5)) 
  
  mutations0all <- length(which( dat$pop$aim3_no_muts==0))
  mutations1all <- length(which( dat$pop$aim3_no_muts==1))
  mutations2all <- length(which( dat$pop$aim3_no_muts==2))
  mutations3all <- length(which(dat$pop$aim3_no_muts==3))
  mutations4all <- length(which(dat$pop$aim3_no_muts==4))
  mutations5all <- length(which(dat$pop$aim3_no_muts==5))

  mutations1plusall <- length(which( dat$pop$aim3_no_muts>=1))
  mutations2plusall <- length(which( dat$pop$aim3_no_muts>=2))
  mutations3plusall <- length(which(dat$pop$aim3_no_muts>=3))
  mutations4plusall <- length(which(dat$pop$aim3_no_muts>=4))

    #coital acts
  
  if(!is.null(dat$discord_coital_df)){
    number_coit_acts <- sum(tapply(dat$discord_coital_df$act_id_couple,
                            dat$discord_coital_df$couple_id,
                            max))
    acts_iev <- length(which(dat$discord_coital_df$iev==1))/2
    percent_iev <-  (acts_iev / number_coit_acts)
    transmission_opps_condom_percent <- (length(which(dat$discord_coital_df$condom==1)) / 
                                   nrow(dat$discord_coital_df) )
     trans_probs_mean <- mean(dat$discord_coital_df$trans_probs)
   
  }else{
    number_coit_acts <- 0
    percent_iev <- NA
    percent_condom <- NA
    trans_probs_mean <- NA
    transmission_opps_condom_percent <- NA
  }
    
  
    #actual calculation of summary stats based on indices and vectors from above
    # and functions for each variable in "popsumm_fxns"
    popsumm_vars=names(dat$popsumm)
  for(ii in 1:length(popsumm_vars)){
    temp_var<-popsumm_vars[ii]
    environment(dat$popsumm_fxns[[ii]])<-environment()
    dat$popsumm[[temp_var]][popsumm_index] <- dat$popsumm_fxns[[ii]]()
   }
  
    #calculation of generic attribute stats
    
    #what percent of alive agents are in each category
    #stat: generic_att_percent_cat_xx (xx=1,..,total number of attributes)
    
    #what percent of alive agents are infected in each category
    #stat: generic_att_percent_inf_cat_xx
    
    #stats for generic attribute values need to be treated separately
    #as the number of attributes may vary between model scenarios
    #note: objects below need to be renamed for clarity 
    temp_length <- length(dat$param$generic_nodal_att_values)
    if(temp_length>1){  
       
       #how many alive agents in each category
       temp_table=table(dat$pop$att1[alive_index])
       #how many alive and infected agents in each category
       temp_table2=table(dat$pop$att1[inf_index])
       #total agents
       sum_temp_table=sum(temp_table)
       #this vector makes sure categories from tables above are
       #arranged in ascending order (necessary if zero agents in a particular 
       #category, which would mean they are missing in tables above
       temp_match=match(names(temp_table),1:temp_length)
       
       
       for(zz in 1:length(temp_match)){
         namevec <- paste("generic_att_percent_cat_",temp_match[zz],sep="")
         dat$popsumm[[namevec]][popsumm_index]=temp_table[zz]/sum_temp_table
       }
       for(zz in 1:temp_length){
         namevec2 <- paste("generic_att_percent_inf_cat_",zz,sep="")
         
         ix1<- which(names(temp_table)==zz)
         ix2<- which(names(temp_table2)==zz)
         if(length(ix2)>0){
           val<- temp_table2[ix2]/temp_table[ix1]  
         }else{val<-0}
         dat$popsumm[[namevec2]][popsumm_index] <- val
       }
    }
  # end of calculating summary stats for generic attributes    
return(dat)
}