[R-sig-eco] Using NADA and rkt for seasonal Mann Kendall tests for censored datasets

[clarkbar]

I have recently put together some code for running trend analyses on censored
datasets with shifting Method Detection Limits (MDLs). I have applied the
NADA package to determine the log mean and log standard deviation; these are
then used to generate log-normal random numbers which then replace the <MDL
observations. From this, I apply the rkt package for seasonal Mann-Kendall
trend test to determine. 
I know the NADA package has cenken, but that has limited covariate and
seasonal capabilities, and I believe my approach below is more conservative. 

I'm sorry if the code is sloppy, but I'm pretty new to all this.. going on
nearly 2 months now. 
I would really appreciate your feedback and ideas. I think my method is
similar to <a
href="http://dugi-doc.udg.edu/bitstream/handle/10256/708/Olea.pdf?sequence=1">Olea's,
but I could never get ILR to cooperate with me. 

Thanks in advance. 

#example data creation - probably not the best way to do this... 
test<-matrix(c(rep(1:7,each=20,1),rep(1:2,each=70), 
         rep(rlnorm(10,-5.6,.2),2), 
         rep(rlnorm(10,-5.6,.3)*1.2,2), 
         rep(rlnorm(10,-5.6,.3)*1.22,2), 
         rep(rlnorm(10,-5.6,.2)*1.25,2), 
         rep(rlnorm(10,-5.6,.3)*1.3,2), 
         rep(rlnorm(10,-5.6,.35)*1.4,2), 
         rep(rlnorm(10,-5.6,.35)*1.4,2)), 
         ncol=3) 
colnames(test)<-c('Year',"Season","Obs") 

test<-as.data.frame(test) 
with(test,rkt(Year,Obs,Season,correct=T,rep='a')) 
with(test,plot(Obs~Year)) 

#function below 
require("NADA") 
require("rkt") 
MDL_fix<-function(value,detection,year,block,CV,month,plot_title) 
  #Observed Value vector, Detection Limit Vector,Block Vector for seasonal
MannKendall,covariate,month,Plot title 
{ 
  YearMonth=year+month/12 
  #using NADA package, determine log-mean of and log-standard deviation of
the observed valued 
  #if(mean(value)<1) {#determine if log or -log should be used (cenfit has
issues with negatives) 
  #lmeanvalue<-{sum(summary(cenmle(value,value<detection))[[9]][1:2])} 
  #lsdvalue<-{summary(cenmle(value,value<detection))[[6]]} 
  #} 
  #if(mean(value)>=1) {#determine if log or -log should be used (cenfit has
issues with negatives) 
  if(length(which(value<detection))/length(value)>.50){ 
    return("ERROR: TOO MANY POINTS BELOW MDL")} 
  lmeanvalue<-{(median(cenfit((log(value)),value<detection))[[1]])} 
    lsdvalue<-{(sd(cenfit(log(value),value<detection)))[[1]]} 
  #} 
  if(length(which(value<detection))/length(value)>.05){ #only pursue MDL
fix if >5% of data censored, otherewise: do 1/2 MDL method 
    #create matrix to fill with Mann-Kendall Output 
    rkt_1<-matrix(1,runs,12) 
    colnames(rkt_1)<- c('2-sided
p-value','Score',"Slope","Var(Score)",'2-sided
p-value_corrected','Var(Score)_corrected', 
                        'Partial Score', 'Partial 2-sided p-value', "Partial
Var(Score)", "Partial 2-sided p-value_corrected" 
                        , 'Partial Var(Score)_corrected',"Tau") 
    for(i in seq(1:runs)) { 
{ 
  xj<-1:length(value) #create random number vector 
  MDL_fixed<-1:length(value) #create final value vector 
  for(j in seq(1:length(value))) { 
{ 
  if(value[j]<detection[j]) #only pursue number generation if the observed
value is below the detection limit 
{ 
    #create a random number (distributed log-normally based on lmean and
lsd) that is below the detection limit 
    repeat { 
      xj[j]<-rlnorm(1,lmeanvalue,lsdvalue) 
      if(xj[j]<detection[j]) {break} 
    } 
  } 
} 
    MDL_fixed[j]<-ifelse(value[j]<detection[j],xj[j],value[j])
#replace values below the detection limit with generated number 
  } 
} 
if(!missing(CV)& !missing(block)) { #if covariate and block present, do
seasonal mann-kendall 
rkt_1[i,]<-as.vector(unlist(rkt(year,MDL_fixed,block,CV,correct=T,rep="a")))
#run MK for each random number generation 
} 
if(missing(block) | missing(CV)){  #if covariate and block mising, do
mann-kendall 
 
rkt_1[i,]<-as.vector(unlist(rkt(year,MDL_fixed,correct=F,rep="a")))
#run MK for each random number generation 
} 
    } 
#plot final iteration of MDL_fixed 
plot(YearMonth,MDL_fixed,
log='y',main=as.character(substitute(plot_title)),xlab="Year",ylab='Concentration(ppm)') 
yaxis<-10^(par('usr')[c(3,4)]) #grab axis size 
#draw first quartile slope and intercept 
abline(a=median(value[which(value>detection)])-fivenum(rkt_1[,"Slope"])[2]*median(year[which(value>detection)]),b=fivenum(rkt_1[,"Slope"])[2],
untf=TRUE, col="red") 
#draw median slope and intercept 
abline(a=median(value[which(value>detection)])-fivenum(rkt_1[,"Slope"])[3]*median(year[which(value>detection)]),b=fivenum(rkt_1[,"Slope"])[3],
untf=TRUE, col="blue") 
#draw third quartile slope and intercept 
abline(a=median(value[which(value>detection)])-fivenum(rkt_1[,"Slope"])[4]*median(year[which(value>detection)]),b=fivenum(rkt_1[,"Slope"])[4],
untf=TRUE, col="green") 
par(new=T) #hold on 
#plot detection limit 
plot(YearMonth,detection, ylim=yaxis,
log='y',type='l',axes=F,xlab='',ylab='',lty=3) 
par(new=F) #hold off 
return(rkt_1) #return matrix of iterations 
  } 
else{ 
  value_half_MDL<-ifelse(value>detection,value,0.5*detection) 
  if(!missing(CV)& !missing(block)) {#if covariate and block present, do
seasonal mann-kendall 
 
rkt_noMDL=as.vector(unlist(rkt(year,value_half_MDL,block,CV,correct=T,rep="a")))#run
MK w/o num geneartion 
  } 
  if(missing(CV)| missing(block)) {#if covariate and block mising, do
mann-kendall 
   
rkt_noMDL=as.vector(unlist(rkt(year,value_half_MDL,correct=F,rep="a")))#run
MK w/o num geneartion 
  } 
  #plot value over time 
  plot(YearMonth,value_half_MDL,
log='y',main=as.character(substitute(plot_title)),xlab="Year",ylab='Concentration(ppm)') 
  yaxis<-10^(par('usr')[c(3,4)])#grab axis size 
  par(new=T) 
  plot(YearMonth,detection, ylim=yaxis,
log='y',type='l',axes=F,xlab='',ylab='',lty=3) 
  par(new=F) 
 
abline(a=median(value[which(value>detection)])-rkt_noMDL[3]*median(year[which(value>detection)]),b=rkt_noMDL[3],
untf=TRUE, col="red") 
  
  return(rkt_noMDL) #return vector of Mann Kendall 
} 
} ##MDL fixing code for log-normal 

runs=100 
summary(with(test,MDL_fix(Obs,detection=c(rep(.005,80),rep(.002,60)),Year,Season,,month=rep(12,140)),"hello"))



--
View this message in context: http://r-sig-ecology.471788.n2.nabble.com/Using-NADA-and-rkt-for-seasonal-Mann-Kendall-tests-for-censored-datasets-tp7578989.html
Sent from the r-sig-ecology mailing list archive at Nabble.com.