R/cc_deaths.R

In corona: Coronavirus ('Rona') Data Exploration

# Generated by LaTeX DogWagger Version 4.0.5 from file <NCTLL_904.tex>
# Date: [2020-9-17 13:16:5] 
# Do NOT edit this file. Edit the LaTeX source!!

# - <Section 5> - 
#' Plot country deaths by week, with various adjustments:
#'
#' Assumes the existence of the data frame stmf containing relevant iso_codes for countries. 
#' The unusual codes GBRTENW and GBR_SCO represent England+Wales and Scotland.
#' You can obtain a list of countries by country_dead('?'), forcing a diagnostic error! 
#' 
#' The columns in the frame stmf are just 'iso_code', 'Year', 'Week', and 'Deaths'. 
#' 
#'   Draws three graphs:
#'   1. Raw data with a linear regression line, over n years; 
#'   2. Data with secular adjustment; 
#'   3. Data adjusted for a 'summer baseline' using the "other n years of data" after
#'      secular adjustment.  
#' 
#' @param pdf default FALSE will not print to PDF
#' @param country Country name 
#' @param save Do we save the data as a CSV
#' @keywords corona deaths 
#' @export 
#' @import ggplot2 
#' @importFrom utils write.csv
#' @importFrom stats lm 
#' @examples 
#' country_dead( 'New Zealand' ) 

# - <Section 6> - 
country_dead <- function ( country='England+Wales', pdf=FALSE, save=FALSE ) 
{ cc <- country_code(country, stmf); 
  mytitle <- paste(country , "Weekly Deaths"); 
  ccname <- paste(cc, '_weekly', sep=''); 
  PdfFilename <- paste(ccname, '_dead.pdf', sep=''); 
  SecularFilename <- paste(ccname, '_secular.pdf', sep=''); 
  AdjustedFilename <- paste(ccname, '_adjusted.pdf', sep=''); 
  MYDEAD <- stmf[ stmf$iso_code==cc, ]; 
  ccRows <- nrow(MYDEAD);

  ## For each row, WE NOW NEED TO insert a Date field based on year+month 
  MinYear <- min(MYDEAD$Year); 
  MaxYear <- max(MYDEAD$Year); 
  MYDEAD$Num <- seq(ccRows);
  startDate <- as.Date(ISOdate(MinYear, 1, 1)); # assumes 1 January [explore, hmm] cf. 1st Monday or Sunday ?
  MYDEAD$Date <- seq.Date(from=startDate, by=7, length.out=ccRows); 
  ccYears <- MaxYear - MinYear; 
if(ccYears < 4)  # pushing our luck, 5 might be better as minimum. 
  { stop( paste("Too few years =", ccYears, "for country", cc) );
  }; 
  YearTop <- ccYears * 52; # assumes first year starts at week 1 [check, hmm]

# - <Section 7> - 
  lowstart <- cntry[ cntry$iso_code==cc, c('lowstart')]; # 'low' is 'summer' 
  lowend <- cntry[ cntry$iso_code==cc, c('lowend')]; # here might check sanity of lowstart, lowend. 
  print( paste('Low range for', country, '=', lowstart, '..', lowend) ); 
if(lowstart < lowend) 
  { summer <- function (M, sStart, sEnd) # N hemisphere
      { M$summer <- as.integer(format(M$Date, '%m')); 
        M$summer <- (M$summer >= sStart) & (M$summer <= sEnd); 
        return(M); 
      }; 
  } else
  { summer <- function (M, sStart, sEnd)
      { M$summer <- as.integer(format(M$Date, '%m')); 
        M$summer <- (M$summer >= sStart) | (M$summer <= sEnd); 
        return(M); 
      }; 
  }; 

# - <Section 8> - 
  secular_fix <- function (F) # F is the frame.
  { mylm <- lm(F$Deaths ~ F$Num);
    yint <- mylm$coefficients[1];
    slp <- mylm$coefficients[2]; 
    WTOT <- nrow(F); 
    F$secular <- F$Deaths * ( 1 + slp*(1+WTOT-F$Num)/yint ); 
    return(F);
  };

  MYDEAD <- secular_fix(MYDEAD); 
  #####################################
  # eye candy with regression line. 
  corona_pdf(PdfFilename, 12, 7, pdf); 
  myplot <- ggplot( MYDEAD, aes(x=.data$Date, y=.data$Deaths)) + 
    labs(title=mytitle, x='Date', y='Deaths') + 
    geom_point() + geom_line() + geom_smooth(method='lm', colour='red', linetype='solid') + 
    scale_x_date(breaks="1 year", date_labels = "%Y") ;
  corona_print(myplot); 
  corona_pdf_off(PdfFilename, pdf, waitline=TRUE);  

  #########################################
  # with secular adjustment:
  corona_pdf(SecularFilename, 12, 7, pdf); 
  secplot <- ggplot( MYDEAD, aes(x=.data$Date, y=.data$secular)) + 
    labs(title=mytitle, x='Date', y='Deaths(adj)') + 
    geom_point() + geom_line() + geom_smooth(method='lm', colour='red', linetype='solid') + 
    scale_x_date(breaks="1 year", date_labels = "%Y") ;
  corona_print(secplot); 
  corona_pdf_off(PdfFilename, pdf, waitline=TRUE);

# - <Section 9> - 
  MYDEAD <- summer(MYDEAD, lowstart, lowend); 
  SumSecAvg <- mean( MYDEAD[ MYDEAD$summer, ]$secular ); 
  MYDEAD$bm <- MYDEAD$Week %% 52;
  weekmean <- tapply(MYDEAD$secular[1:YearTop], MYDEAD$bm[1:YearTop], mean);
  MYDEAD$expected <- MYDEAD$adjust <- weekmean[ 1+MYDEAD$bm ];
  qititle <- paste('Seasonally Adjusted Deaths:', country); 

  # next, remove own value from mean: 
  MYDEAD$adjust[1:YearTop] <- (ccYears*MYDEAD$adjust[1:YearTop] 
                               - MYDEAD$secular[1:YearTop])/(ccYears-1); 
  ## ^ do not adjust final values over YearTop i.e. in 2020, as they haven't been included! 
  # Now normalise to SumSecAvg, the summer secular average. 
  MYDEAD$adjust = MYDEAD$adjust/SumSecAvg; 
  MYDEAD$seasonal <- MYDEAD$secular/MYDEAD$adjust; # finally, adjust all. 

  corona_pdf(AdjustedFilename, 12, 7, pdf); 
  mychart <- qic(x=MYDEAD$Date, y=MYDEAD$seasonal, chart='i',
          title=qititle, xlab='Year', ylab='Deaths'); 
  # print(mychart); 
  corona_print(mychart, islegend=FALSE);
  corona_pdf_off(AdjustedFilename, pdf); 
if(save)
  { CsvFilename <- paste(ccname, 'deaths.csv', sep=''); 
    write.csv(MYDEAD, file=CsvFilename);
  }; 
}
# -END OF FILE-