Function_SNI_Kelly2011.R

CharSNI = function(CharData,BandWidth) {
  ## Calculate signal-to-noise index (SNI) for a charcoal record (see Kelly et al., 2011).
  #
  # CharData = Matrix of input data with one row per sample, containing:
  #   Column 1: age associated with the sample (yr)
  #   Column 2: charcoal accumulation rate (CHAR) of the sample (pieces/cm^2/yr)
  #   Column 3: threshold value (pieces/cm^2/yr)
  #
  # BandWidth = Width of moving window for computing SNI 
  #
  #
  # This function computes SNI as described in Kelly et al. 2010.  Note that
  # your data must be interpolated to constant sample resolution (yr/sample)
  # before input to the function. The function makes no assumption about
  # prior analysis on the input CHAR series, i.e. any background and
  # threshold methods may be used.  However, input data should still be
  # consistent with the interpretation that a CHAR value (column 2) greater
  # than the corresponding threshold value (column 3) is a "signal" sample,
  # whereas a CHAR value below the threshold is "noise".  Refer to Kelly et
  # al. 2010 for details and discussion.
  #
  # SNI_output is a data list containing the computed SNI and related
  # data, with one row for each row in the input variable CharData:
  #   
  #   $SNI    = the SNI computed for each sample
  #   $winInd = indexes of the first and last samples included in each moving 
  #     window.  E.g. SNI_output$winInd(X) == [A, B] indicates that the 
  #     moving window used to calculate SNI for the Xth sample contained all 
  #     samples between A and B, inclusive.
  #   $popN   = the CHAR values of all samples in the noise (N) population 
  #     (samples in the moving window with CHAR below threshold)
  #   $popS   = the CHAR values of all samples in the signal (S) population 
  #     (samples in the moving window with CHAR at or above threshold)
  #   $meanN  = mean CHAR of the samples in popN
  #   $stdN   = standard deviation of the samples in popN
  #   $CF     = the "correction factor" used in computing SNI.  Equal to
  #     (v - 2)/v ,where v is the number of samples in popN
  
  # load required library
  library(stats)
  
  # Data setup
  ages = CharData[,1];
  CHAR = CharData[,2];
  thresh = CharData[,3];
  
  if( sd(diff(ages))==0 ) { # If resolution (yr/sample) is constant
    r = mean(diff(ages)); # r = constant sampling resolution
  } else {
    print('Fatal error: Sampling resolution (yr/sample) must be constant.');
    return();
  }
  
  # Preallocate some space   
  SNI_output 				= list();
  SNI_output$SNI 		= NA*numeric(length(ages));
  SNI_output$winInd	= matrix(data=NA,nrow=length(ages),ncol=2);
  SNI_output$popN 	= vector(mode="list",length=length(ages));
  SNI_output$popS 	= vector(mode="list",length=length(ages));
  SNI_output$meanN 	= NA*numeric(length(ages));
  SNI_output$stdN 	= NA*numeric(length(ages));
  SNI_output$CF 		= NA*numeric(length(ages));
  rawSNI 						= NA*numeric(length(ages));
  
  for(i in 1:length(ages)) { # Perform calculations at each sample age
    # Calculate window indexes
    if( i < round(0.5*(BandWidth/r))+1 ) {  # Samples near beginning (moving window truncated)
      SNI_output$winInd[i,] = c(1, ( i + round(0.5*(BandWidth/r)) )); 
    } else {
      if( i > length(ages)-round(0.5*(BandWidth/r)) ) { # Samples near end
        SNI_output$winInd[i,] = c( (i - round(0.5*(BandWidth/r))), length(ages) );  
      } else { 
        SNI_output$winInd[i,] = c( (i-round(0.5*(BandWidth/r))), (i+round(0.5*(BandWidth/r))) );
      }
    }
    
    # In each moving window...  
    
    # Get CHAR for samples in window (X is entire window population)
    X = CHAR[SNI_output$winInd[i,1]:SNI_output$winInd[i,2]]; 
    
    # inS, inN: boolean indexes of samples counted as S & N, respectively
    inS = (X >= thresh[SNI_output$winInd[i,1]:SNI_output$winInd[i,2]] );
    inN = !inS;
    
    # Fill in S & N populations, means, stds
    SNI_output$popS[[i]] = X[inS];
    SNI_output$popN[[i]] = X[inN];  
    SNI_output$meanN[i] = mean(X[inN]);
    SNI_output$stdN[i] = sd(X[inN]);
    
    # Calculate correction factor
    v = length(SNI_output$popN[[i]]); # d.f. of N
    SNI_output$CF[i] = (v - 2) / v; # SNI = Z * (v-2)/v
    
    # Calculate raw SNI (unsmoothed)
    if( length(SNI_output$popS[[i]])>0 ) {
      rawSNI[i] = (mean( (SNI_output$popS[[i]] - SNI_output$meanN[i]) ) 
                   / SNI_output$stdN[i]) * SNI_output$CF[i];
    } else {
      rawSNI[i] = 0; # SNI = 0 by definition when no samples exceed threshold
    }
    
  }
  
  # Smooth raw values to obtain final SNI
  SNI_output$SNI = lowess(ages,rawSNI,f=(BandWidth/r)/length(ages),iter=0)$y;
  
  return(SNI_output);
  
}