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.
Using NADA and rkt for seasonal Mann Kendall tests for censored datasets
2 messages · clarkbar, Philippi, Tom
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-ecology/attachments/20140724/788a1d12/attachment.pl>