################################################################################
#
# The 2022 Europe report of the Lancet Countdown on health and climate change: 
# towards a climate resilient future
# https://doi.org/10.1016/S2468-2667(22)00197-9
#
# Indicator 1.1.4: Heat Related Mortality
# (C) Joan Ballester (joan.ballester@isglobal.org), Hicham Achebak
# March 2022
# 
################################################################################

rm( list = ls() );
cat("\014");

# Required Libraries
suppressMessages( library(data.table) ); # Required Functions: fread
suppressMessages( library(dlnm) ); # Required Functions: crossbasis
suppressMessages( library(splines) ); # Required Functions: ns
suppressMessages( library(mixmeta) ); # Required Functions: mixmeta
suppressMessages( library(tsModel) ); # Required Functions: Lag
suppressMessages( library(MASS) ); # Required Functions: mvrnorm



################################################################################
# Definition of the General Parameters
################################################################################

print("= Definition of the General Parameters =");

# Period of Time
YEAR1 = 2015;
YEAR2 = 2020;

# Degrees of Freedom per Year for the Seasonal and Long-Term Trends
DF_SEAS = 8;

# Exposure-Response: Percentiles of the Temperature Knots
VAR_PRC = c(10,75,90) / 100;

# Exposure-Response: Type of the Spline
VAR_FUN = "ns";

# Lag-Response: Minimum and Maximum Lags (in Weeks)
MIN_LAG = 0;
MAX_LAG = 3;

# Cumulative Exposure-Response: Percentiles for the Predictions
PRED_PRC = c( seq(  0.0,   1.0, 0.1 ),
              seq(  1.5,   5.0, 0.5 ),
              seq(  6.0,  94.0, 1.0 ),
              seq( 95.0,  98.5, 0.5 ),
              seq( 99.0, 100.0, 0.1 ) ) / 100;

# Lower and Upper Percentiles for the Calculation of the Minimum Mortality Temperature
MIN_PMMT =  5 / 100;
MAX_PMMT = 99 / 100;



################################################################################
# Preparation of the Data
################################################################################

print("= Preparation of the Data =");

# Reading the Datatable
DATATABLE = fread( "data.csv" );

# Reading the Metatable
METATABLE = fread( "metadata.csv", encoding = "Latin-1" );

# Creating the Vector of Regional NUTS Codes
vREG = METATABLE$location;
nREG = length(vREG);

# Creating the Vector of Regional Names
vREG_NAME = METATABLE$name;
rm(METATABLE);

# Creating the Datalists
DATALIST = lapply( vREG, function(x) DATATABLE[ DATATABLE$location == x, ] );
names(DATALIST) = vREG;
rm(DATATABLE);

# Adding the Variable Week-of-Period
for( iREG in 1:nREG ){
  DATALIST[[iREG]]$wop = 1:length(DATALIST[[iREG]]$woy);
}
rm(iREG);



################################################################################
# Calculation of the Location-Specific Associations
################################################################################

print("= Calculation of the Location-Specific Associations =");

# Formula of the Model
#   length(wop)              :: Number of Weeks of the Period
#   length(wop) * 7          :: Number of Days of the Period
#   length(wop) * 7 / 365.25 :: Number of Years of the Period
FORMULA = mort ~ ns( wop, df = round( DF_SEAS * length(wop) * 7 / 365.25 ) ) + CROSS_BASIS;

# Reduced Coefficients
COEF_MODEL = matrix( data = NA, nREG, length(VAR_PRC) + 1, dimnames = list(vREG) );

# Reduced Covariance Matrices
VCOV_MODEL = vector( "list", nREG );
names(VCOV_MODEL) = vREG;

# Cumulative Exposure-Response Before the Meta-Analysis
CROSS_PRED_NOBLUP = vector( "list", nREG );
names(CROSS_PRED_NOBLUP) = vREG;

for( iREG in 1:nREG ){
  print( paste0( "   Region ", iREG, ": ", vREG_NAME[iREG], " (", vREG[iREG], ")" ) );
  
  # Cross-Basis of Temperature
  CROSS_BASIS = crossbasis( DATALIST[[iREG]]$temp,
                            c( MIN_LAG, MAX_LAG ),
                            argvar = list( fun = VAR_FUN,
                                           knots = quantile( DATALIST[[iREG]]$temp, VAR_PRC, na.rm = TRUE ),
                                           Boundary.knots = range( DATALIST[[iREG]]$temp, na.rm = TRUE ) ),
                            arglag = list( fun = "integer" ) );

  # Fitting of the Model
  GLM_MODEL = glm( formula = FORMULA, DATALIST[[iREG]], family = quasipoisson, na.action = "na.exclude" );
  
  # Cumulative Exposure-Response without Specifying the Centring Temperature
  suppressMessages( CROSS_PRED_NOBLUP[[iREG]] <- crosspred( CROSS_BASIS, GLM_MODEL, at = quantile( DATALIST[[iREG]]$temp, PRED_PRC, na.rm = TRUE ) ) );
  
  # Vector of Local Minima of the Cumulative Exposure-Response
  iMMT = 1 + which( diff( sign( diff( CROSS_PRED_NOBLUP[[iREG]]$allRRfit ) ) ) == 2 );
  
  # Minimum Mortality Temperature
  if( length(iMMT) > 0 ){
    # If There Is at Least One Local Minimum, the Minimum Mortality Temperature Is the Local Minimum with the Lowest Relative Risk
    MMT = CROSS_PRED_NOBLUP[[iREG]]$predvar[ iMMT[ which.min( CROSS_PRED_NOBLUP[[iREG]]$allRRfit[iMMT] ) ] ];
  }else{
    # If There Are No Local Minima, the Minimum Mortality Temperature Is the Lowest Relative Risk within the Predefined Temperature Percentile Range [MIN_PMMT,MAX_PMMT]
    MMT = CROSS_PRED_NOBLUP[[iREG]]$predvar[ which( PRED_PRC == MIN_PMMT ) - 1 + which.min( CROSS_PRED_NOBLUP[[iREG]]$allRRfit[ which( PRED_PRC == MIN_PMMT ) : which( PRED_PRC == MAX_PMMT ) ] ) ];
  }
  rm(iMMT);
  
  # Cumulative Exposure-Response Specifying the Centring Temperature
  CROSS_PRED_NOBLUP[[iREG]] = crosspred( CROSS_BASIS, GLM_MODEL, at = quantile( DATALIST[[iREG]]$temp, PRED_PRC, na.rm = TRUE ), bylag = 1, cen = MMT );
  
  # Reduced Coefficients and Covariance Matrices
  REDUCED = crossreduce( CROSS_BASIS, GLM_MODEL, cen = MMT );
  COEF_MODEL[iREG,] = coef( REDUCED );
  VCOV_MODEL[[iREG]] = vcov( REDUCED );
  rm(REDUCED);
  
  rm(CROSS_BASIS,GLM_MODEL, MMT);
  
}
rm(FORMULA,iREG);



################################################################################
# Calculation of the Meta-Analysis and the Best Linear Unbiased Predictions
################################################################################

print("= Calculation of the Meta-Analysis and the Best Linear Unbiased Predictions =");

# Minimum Mortality Temperature
MMT = array( NA, dim = nREG, dimnames = list(vREG) );

# Cumulative Exposure-Response After the Meta-Analysis
CROSS_PRED_BLUP = vector( "list", nREG );
names(CROSS_PRED_BLUP) = vREG;

# Multivariate Meta-Analysis of the Reduced Coefficients
MULTIVAR = mixmeta( COEF_MODEL ~ 1, VCOV_MODEL, data = data.frame( vREG = vREG ), control = list( showiter = FALSE ) );

# Best Linear Unbiased Predictions
BLUP = blup( MULTIVAR, vcov = TRUE );

for( iREG in 1:nREG ){
  
  # Basis of Temperature
  BASIS_VAR = onebasis( quantile( DATALIST[[iREG]]$temp, PRED_PRC, na.rm = TRUE ),
                        fun = VAR_FUN,
                        knots = quantile( DATALIST[[iREG]]$temp, VAR_PRC, na.rm = TRUE ),
                        Boundary.knots = range( DATALIST[[iREG]]$temp, na.rm = TRUE ) );
  
  # Vector of Percentile Temperatures for the Predictions
  PRED_VAR = quantile( DATALIST[[iREG]]$temp, PRED_PRC, na.rm = TRUE );
  
  # Cumulative Exposure-Response without Specifying the Centring Temperature
  suppressMessages( PRED_MORT <- crosspred( BASIS_VAR, coef = BLUP[[iREG]]$blup, vcov = BLUP[[iREG]]$vcov, model.link = "log", at = PRED_VAR ) );
  
  # Vector of Local Minima of the Cumulative Exposure-Response
  iMMT = 1 + which( diff( sign( diff( PRED_MORT$allRRfit ) ) ) == 2 );
  
  # Minimum Mortality Temperature
  if( length(iMMT) > 0 ){
    # If There Is at Least One Local Minimum, the Minimum Mortality Temperature Is the Local Minimum with the Lowest Relative Risk
    MMT[iREG] = PRED_VAR[ iMMT[ which.min( PRED_MORT$allRRfit[iMMT] ) ] ];
  }else{
    # If There Are No Local Minima, the Minimum Mortality Temperature Is the Lowest Relative Risk within the Predefined Temperature Percentile Range [MIN_PMMT,MAX_PMMT]
    MMT[iREG] = PRED_VAR[ which( PRED_PRC == MIN_PMMT ) - 1 + which.min( PRED_MORT$allRRfit[ which( PRED_PRC == MIN_PMMT ) : which( PRED_PRC == MAX_PMMT ) ] ) ];
  }
  rm(PRED_VAR,PRED_MORT,iMMT);
  
  # Cumulative Exposure-Response Specifying the Centring Temperature
  CROSS_PRED_BLUP[[iREG]] = crosspred( BASIS_VAR, coef = BLUP[[iREG]]$blup, vcov = BLUP[[iREG]]$vcov, model.link = "log", at = quantile( DATALIST[[iREG]]$temp, PRED_PRC, na.rm = TRUE ), cen = MMT[iREG] );
  rm(BASIS_VAR);
  
}
rm(COEF_MODEL,VCOV_MODEL, MULTIVAR, iREG);



################################################################################
# Calculation of the Attributable Mortality
################################################################################

print("= Calculation of the Attributable Mortality =");

# Vector of Time Periods: Overall Period and Individual Years
vPER = c( paste0( YEAR1 , "-", YEAR2 ), YEAR1:YEAR2 );
nPER = length(vPER);

# Vector of Temperature Ranges
vRNG = c( "Total", "Total Cold", "Total Heat", "Moderate Cold", "Moderate Heat", "Extreme Cold", "Extreme Heat" );
nRNG = length(vRNG);

# Number of Simulations
nSIM = 10^3;

# Values and Simulations of Attributable Number for Individual Regions
ATT_VAL = array( NA, dim = c( nREG, nPER, nRNG       ), dimnames = list( vREG, vPER, vRNG         ) );
ATT_SIM = array( NA, dim = c( nREG, nPER, nRNG, nSIM ), dimnames = list( vREG, vPER, vRNG, 1:nSIM ) );

# Attributable Number (AN, in Deaths), Fraction (AF, in %) and Incidence (AI, in Daily Deaths per Million) for Individual Regions (REG) and the Whole Domain (TOT)
AN_REG = AF_REG = AI_REG = array( NA, dim = c( nREG, nPER, nRNG, 3 ), dimnames = list(    vREG  , vPER, vRNG, c( "att_val", "att_low", "att_upp" ) ) );
AN_TOT = AF_TOT = AI_TOT = array( NA, dim = c(   1 , nPER, nRNG, 3 ), dimnames = list( "Germany", vPER, vRNG, c( "att_val", "att_low", "att_upp" ) ) );

# Mortality (MORT) and Population (POPU) for Individual Regions (REG) and the Whole Domain (TOT)
MORT_REG = POPU_REG = array( NA, dim = c( nREG, nPER ), dimnames = list(    vREG  , vPER ) );
MORT_TOT = POPU_TOT = array( NA, dim = c(   1 , nPER ), dimnames = list( "Germany", vPER ) );

# Regional Thresholds for the Temperature Ranges
TMIN = P025 = P975 = TMAX = array( NA, dim = nREG, dimnames = list(vREG) );
for( iREG in 1:nREG ){
  TMIN[iREG] =      min( DATALIST[[iREG]]$temp,        na.rm = TRUE );
  P025[iREG] = quantile( DATALIST[[iREG]]$temp, 0.025, na.rm = TRUE );
  P975[iREG] = quantile( DATALIST[[iREG]]$temp, 0.975, na.rm = TRUE );
  TMAX[iREG] =      max( DATALIST[[iREG]]$temp,        na.rm = TRUE );
}
rm(iREG);

# 1/2: Loop for Individual Regions
for( iREG in 1:nREG ){
  
  # Basis of Temperature
  BASIS_VAR = onebasis( DATALIST[[iREG]]$temp,
                        fun = VAR_FUN,
                        knots = quantile( DATALIST[[iREG]]$temp, VAR_PRC, na.rm = TRUE ),
                        Boundary.knots = range( DATALIST[[iREG]]$temp, na.rm = TRUE ) );

  # Minimum Mortality Temperature for the Centring of the Basis of Temperature
  BASIS_MMT = onebasis( MMT[iREG],
                        fun = VAR_FUN,
                        knots = quantile( DATALIST[[iREG]]$temp, VAR_PRC, na.rm = TRUE ),
                        Boundary.knots = range( DATALIST[[iREG]]$temp, na.rm = TRUE ) );

  # Basis of Temperature After Centring at the Minimum Mortality Temperature
  BASIS_CEN = scale( BASIS_VAR, center = BASIS_MMT, scale = FALSE );
  rm(BASIS_VAR,BASIS_MMT);
  
  # Matrix of Deaths with Lags
  LAGGED_MORT_MATRIX = Lag( DATALIST[[iREG]]$mort, -MIN_LAG:-MAX_LAG );
  
  # Average of Deaths Across Lags
  LAGGED_MORT_VECTOR = rowMeans( LAGGED_MORT_MATRIX );
  
  # Attributable Number: Weekly Time Series
  ATT_NUM_TS_REF = ( 1 - exp( -BASIS_CEN %*% BLUP[[iREG]]$blup ) ) * LAGGED_MORT_VECTOR;
  
  # Perturbed Coefficients of the BLUPs (Dimensions: Coefficients x Simulations)
  set.seed( iREG );
  COEF_SIM = t( mvrnorm( nSIM, BLUP[[iREG]]$blup, BLUP[[iREG]]$vcov ) );
  
  # Attributable Number: Matrix of Perturbed Weekly Time Series (Dimensions: Time x Simulations)
  ATT_NUM_TS_SIM = ( 1 - exp( -BASIS_CEN %*% COEF_SIM ) ) * LAGGED_MORT_VECTOR;
  rm(BASIS_CEN,COEF_SIM);
  
  # Loop for the Vector of Time Periods
  for( iPER in 1:nPER ){
    
    # Selection of Time Period Data
    if( iPER == 1 ){ vTIM = which( YEAR1 <= DATALIST[[iREG]]$year & DATALIST[[iREG]]$year <= YEAR2 ); } # All Years
    else           { vTIM = which(               DATALIST[[iREG]]$year == vPER[iPER]               ); } # Individual Years
    if( length(vTIM) <= 0 ){ stop("ERROR: Invalid Time Period Data !!!"); }
    
    # Mortality and Population of the Region and Time Period
    MORT_REG[iREG,iPER] =  sum( DATALIST[[iREG]]$mort[vTIM], na.rm = TRUE );
    POPU_REG[iREG,iPER] = mean( DATALIST[[iREG]]$popu[vTIM], na.rm = TRUE );
    
    # Loop for the Vector of Temperature Ranges
    for( iRNG in 1:nRNG ){
      
      # Indicator of Temperature Range      
      if     ( iRNG == 1 ){ vRNG_THRES =                                           1:length(vTIM)                                    ; } # Total
      else if( iRNG == 2 ){ vRNG_THRES = which(                                            DATALIST[[iREG]]$temp[vTIM] <  MMT[iREG] ); } # Total Cold
      else if( iRNG == 3 ){ vRNG_THRES = which(  MMT[iREG] < DATALIST[[iREG]]$temp[vTIM]                                            ); } # Total Heat
      else if( iRNG == 4 ){ vRNG_THRES = which( P025[iREG] < DATALIST[[iREG]]$temp[vTIM] & DATALIST[[iREG]]$temp[vTIM] <  MMT[iREG] ); } # Moderate Cold
      else if( iRNG == 5 ){ vRNG_THRES = which(  MMT[iREG] < DATALIST[[iREG]]$temp[vTIM] & DATALIST[[iREG]]$temp[vTIM] < P975[iREG] ); } # Moderate Heat
      else if( iRNG == 6 ){ vRNG_THRES = which(                                            DATALIST[[iREG]]$temp[vTIM] < P025[iREG] ); } # Extreme Cold
      else if( iRNG == 7 ){ vRNG_THRES = which( P975[iREG] < DATALIST[[iREG]]$temp[vTIM]                                            ); } # Extreme Heat
      else                { stop("ERROR: Invalid Temperature Range !!!"); }
      if( length(vRNG_THRES) > 0 ){
        
        # Attributable Number: Value by Temperature Range
        ATT_VAL[iREG,iPER,iRNG ] =     sum( ATT_NUM_TS_REF[ vTIM[vRNG_THRES]               ], na.rm = TRUE ) * sum( rowMeans( LAGGED_MORT_MATRIX[ vTIM, , drop=FALSE ], na.rm = TRUE ), na.rm = TRUE ) /
                                                                                                               sum(           LAGGED_MORT_VECTOR[ vTIM               ]                , na.rm = TRUE );
        
        # Attributable Number: Simulations by Temperature Range
        ATT_SIM[iREG,iPER,iRNG,] = colSums( ATT_NUM_TS_SIM[ vTIM[vRNG_THRES], , drop=FALSE ], na.rm = TRUE ) * sum( rowMeans( LAGGED_MORT_MATRIX[ vTIM, , drop=FALSE ], na.rm = TRUE ), na.rm = TRUE ) /
                                                                                                               sum(           LAGGED_MORT_VECTOR[ vTIM               ]                , na.rm = TRUE );
        
        # Attributable Number: Value and Confidence Interval
        AN_REG[iREG,iPER,iRNG, 1 ] =           ATT_VAL[iREG,iPER,iRNG ]                  ;
        AN_REG[iREG,iPER,iRNG,2:3] = quantile( ATT_SIM[iREG,iPER,iRNG,], c(0.025,0.975) );
        
        # Attributable Fraction
        AF_REG[iREG,iPER,iRNG,] = 10^2 * AN_REG[iREG,iPER,iRNG,] / MORT_REG[iREG,iPER];
        
        # Attributable Incidence (in Daily Deaths per Million)
        #   length(vTIM)     :: Number of Weeks Included in the Calculation
        #   length(vTIM) * 7 :: Number of Days Included in the Calculation
        AI_REG[iREG,iPER,iRNG,] = 10^6 * AN_REG[iREG,iPER,iRNG,] / ( POPU_REG[iREG,iPER] * length(vTIM) * 7 );
        
      }else{
        ATT_VAL[iREG,iPER,iRNG ] = 0;
        ATT_SIM[iREG,iPER,iRNG,] = 0;
        AN_REG[iREG,iPER,iRNG,] = 0;
        AF_REG[iREG,iPER,iRNG,] = 0;
        AI_REG[iREG,iPER,iRNG,] = 0;
      }
      rm(vRNG_THRES);
      
    }
    rm(vTIM, iRNG);
    
  }
  rm(LAGGED_MORT_MATRIX,LAGGED_MORT_VECTOR, ATT_NUM_TS_REF,ATT_NUM_TS_SIM, iPER);
  
}
rm(MMT,BLUP, TMIN,P025,P975,TMAX, iREG)

# 2/2: Whole Domain

# Mortality and Population of the Whole Domain and Time Period
MORT_TOT[1,] = colSums( MORT_REG );
POPU_TOT[1,] = colSums( POPU_REG );

# Attributable Number: Value and Confidence Interval
AN_TOT[1,,,1] =        apply( ATT_VAL, c(2,3  ), sum )                           ;
AN_TOT[1,,,2] = apply( apply( ATT_SIM, c(2,3,4), sum ), c(1,2), quantile, 0.025 );
AN_TOT[1,,,3] = apply( apply( ATT_SIM, c(2,3,4), sum ), c(1,2), quantile, 0.975 );

# Loop for the Vector of Time Periods
for( iPER in 1:nPER ){
  
  # Selection of Time Period Data
  if( iPER == 1 ){ vTIM = which( YEAR1 <= DATALIST[[1]]$year & DATALIST[[1]]$year <= YEAR2 ); } # All Years
  else           { vTIM = which(              DATALIST[[1]]$year == vPER[iPER]             ); } # Individual Years
  if( length(vTIM) <= 0 ){ stop("ERROR: Invalid Time Period Data !!!"); }
  
  # Attributable Fraction
  AF_TOT[1,iPER,,] = 10^2 * AN_TOT[1,iPER,,] / MORT_TOT[1,iPER];
  
  # Attributable Incidence (in Daily Deaths per Million)
  #   length(vTIM)     :: Number of Weeks Included in the Calculation
  #   length(vTIM) * 7 :: Number of Days Included in the Calculation
  AI_TOT[1,iPER,,] = 10^6 * AN_TOT[1,iPER,,] / ( POPU_TOT[1,iPER] * length(vTIM) * 7 );
  rm(vTIM);
  
}
rm(iPER);

rm(ATT_VAL,ATT_SIM, MORT_REG,POPU_REG, MORT_TOT,POPU_TOT);



################################################################################
# Cumulative ERF of Individual Regions Before and After the Meta-Analysis
################################################################################

print("= Cumulative Exposure-Response Function of Individual Regions =");

# Before the Meta-Analysis

pdf( "PLOT_CUMU_ERF_NOBLUP.pdf", width = 15, height = 15 );
layout( matrix( seq(4*4), nrow = 4, byrow = TRUE ) );
par( mex = 0.8, mgp = c(2.5,1,0), las = 0 );

for( iREG in 1:nREG ){
  plot   (    CROSS_PRED_NOBLUP[[iREG]]$predvar                                            ,    CROSS_PRED_NOBLUP[[iREG]]$allRRfit                                              , col = rgb(0.00,0.00,0.00    ), lwd = 4, lty = 1, type = "l", xlab = expression( paste( "Temperature (", degree, "C)" ) ), ylab = "Relative Risk", main = paste0( vREG_NAME[iREG], " (", vREG[iREG], ")" ), ylim = c( 1, max( CROSS_PRED_NOBLUP[[iREG]]$allRRfit ) ), axes = T );
  polygon( c( CROSS_PRED_NOBLUP[[iREG]]$predvar, rev( CROSS_PRED_NOBLUP[[iREG]]$predvar ) ), c( CROSS_PRED_NOBLUP[[iREG]]$allRRlow, rev( CROSS_PRED_NOBLUP[[iREG]]$allRRhigh ) ), col = rgb(0.00,0.00,0.00,1/3), border = FALSE );
  lines  (    CROSS_PRED_NOBLUP[[iREG]]$predvar                                            ,    CROSS_PRED_NOBLUP[[iREG]]$allRRfit                                              , col = rgb(0.00,0.00,0.00    ), lwd = 4, lty = 1 );
  abline( h = 1, lwd = 1, lty = 1, col = rgb(0.00,0.00,0.00) );
}
rm(iREG);

dev.off();

# After the Meta-Analysis

pdf( "PLOT_CUMU_ERF_BLUP.pdf", width = 15, height = 15 );
layout( matrix( seq(4*4), nrow = 4, byrow = TRUE ) );
par( mex = 0.8, mgp = c(2.5,1,0), las = 0 );

for( iREG in 1:nREG ){
  plot   (    CROSS_PRED_BLUP[[iREG]]$predvar                                          ,    CROSS_PRED_BLUP[[iREG]]$allRRfit                                            , col = rgb(0.00,0.00,0.00    ), lwd = 4, lty = 1, type = "l", xlab = expression( paste( "Temperature (", degree, "C)" ) ), ylab = "Relative Risk", main = paste0( vREG_NAME[iREG], " (", vREG[iREG], ")" ), ylim = c( 1, max( CROSS_PRED_BLUP[[iREG]]$allRRfit ) ), axes = T );
  polygon( c( CROSS_PRED_BLUP[[iREG]]$predvar, rev( CROSS_PRED_BLUP[[iREG]]$predvar ) ), c( CROSS_PRED_BLUP[[iREG]]$allRRlow, rev( CROSS_PRED_BLUP[[iREG]]$allRRhigh ) ), col = rgb(0.00,0.00,0.00,1/3), border = FALSE );
  lines  (    CROSS_PRED_BLUP[[iREG]]$predvar                                          ,    CROSS_PRED_BLUP[[iREG]]$allRRfit                                            , col = rgb(0.00,0.00,0.00    ), lwd = 4, lty = 1 );
  abline( h = 1, lwd = 1, lty = 1, col = rgb(0.00,0.00,0.00) );
}
rm(iREG);

dev.off()

#
