######################################################################################################### swaveqexMerge <- function(cdatin,qwstnum,ddstnum,yrbeg,yrend,getdd="WD",minss=3,runname="",redconc=0) { # # merges pesticide data with daily discharge data and prepares # input file for swaveqexFit # # cdatin is a data frame with the pesticide concentration data # - the first column should be the station number (character) # - the second column should be the date, in "yyyy-mm-dd" format (character) # - the concentration value should be a column named "final_value" (numeric) # - the remark should be in a column names "final_remark" (character) # # qwstnum is the station number to analyze # ddstnum is the station number for daily discharge (usually the same as qwstnum) # yrbeg and yrend are the beginning and ending calendar years for analysis (numeric) # # getdd specifies where to get the daily discharge data # getdd="WD" (the default) uses the waterData package to get streamflow # data from the USGS web server # getdd="File" reads daily discharge from the tab-delimited text file # called "dd_ddstnum.txt" in the default directory # minss is the minimum spacing between samples, in days. Samples are thinned, if necessary, # so that there is at least minss days between samples. minss must be an integer # between 1 and 15 if(minss<1) {minss <- 1 } else if (minss>15){minss <- 15 } else {minss <- floor(minss)} # runname is an optional character string for naming the output files produced by # running swaveqexFit. For example, if runname="_Fred_", then the plot output file # will be called "Plots_Fred_XXXXX.pdf" and the conditional simulation file will be # called "CSIM_Fred_XXXXX.txt" where XXXXX=qwstnum # # redconc ("reduce concentration") is a non-negative constant to subtract from concentration before # fitting the model. redconc should be <= (minimum observed concentration -0.001) # # Requires library waterData # # select pesticide data for specified site and # download daily discharge data using importDVs from the waterData library cdatxx <- cdatin[cdatin[,1]==qwstnum,] cdatxx <- data.frame(cdatxx[,2],cdatxx[,"final_remark"],cdatxx[,"final_value"]) names(cdatxx) <- c("date","crem","cval") #### Subtract constant from observed concentrations (if redconc is not zero) #### First make sure that redconc is less than the (minimum concentration - 0.001) if(redconc> (min(cdatxx[,3]) - 0.001)) {redconc <- min(cdatxx[,3])-0.001 } cdatxx[,3] <- cdatxx[,3]-redconc yrxx <- as.numeric(substring(as.character(cdatxx[,1]),1,4)) if(min(yrxx)>yrbeg) {yrbeg <- min(yrxx) } if(max(yrxx)=yrbeg & yr<=yrend,] qdate <- as.character(qin[,1]) qval <- qin[,2] qval[qval<1] <- 1 yr <- as.numeric(substring(qdate,1,4)) mo <- as.numeric(substring(qdate,6,7)) da <- as.numeric(substring(qdate,9,10)) qday <- rep(NA,length(qdate)) cumd <- c(0,31,59,90,120,151,181,212,243,273,304,334) mday <- c(31,28,31,30,31,30,31,31,30,31,30,31) chmo <- c("01","02","03","04","05","06","07","08","09","10","11","12") chday <- c(paste("0",as.character(1:9),sep=""),as.character(10:31)) for (j in 1:12) { qday[mo==j] <- (yr[mo==j]-yrbeg)*365 + cumd[j]+da[mo==j] } qdaily <- rep(NA,365*(yrend-yrbeg+1)) qdaily[qday] <- qval ntot <- 365*(yrend-yrbeg+1) qdate <- rep(" ",ntot) for (i in yrbeg:yrend) { for (j in 1:12) { for (k in 1:mday[j]) { qdate[(i-yrbeg)*365+cumd[j]+k] <- paste(as.character(i),"-",chmo[j],"-",chday[k],sep="") } } } # compute flow anomalies qmean <- mean(log10(qdaily),na.rm=T) qmc <- log10(qdaily)- qmean qanom30 <- rep(0,nall) for (i in 30:nall) {qanom30[i] <- mean(qmc[(i-29):i])} qanom30[1:29] <- qanom30[30] qanom1 <- qmc-qanom30 datxxx <- data.frame(qdate,1:nall,rep(1:365,rlen),qanom30,qanom1) datxxx[,1] <- as.character(datxxx[,1]) names(datxxx) <- c('date','jday','sday','qanom30','qanom1') # # prepare concentration data # chdate2 <- as.character(cdatxx[,1]) yr2 <- as.numeric(substring(chdate2,1,4)) mo2 <- as.numeric(substring(chdate2,6,7)) da2 <- as.numeric(substring(chdate2,9,10)) leapdays <- (mo2==2 & da2==29) cdatxx <- cdatxx[!leapdays & yr2>=yrbeg & yr2<=yrend,] chdate2 <- as.character(cdatxx[,1]) yr2 <- as.numeric(substring(chdate2,1,4)) mo2 <- as.numeric(substring(chdate2,6,7)) da2 <- as.numeric(substring(chdate2,9,10)) jday2 <- (yr2-yrbeg)*365+ cumd[mo2] + da2 sday2 <- cumd[mo2] + da2 #### thin to at least minss days between samples nx2 <- length(jday2) pkthin <- rep(F,nx2) pkthin[1] <- T jdaycum <- jday2[1] for (ii in 2:nx2) { if ((jday2[ii]-jdaycum)>(minss-1)) { pkthin[ii] <- T jdaycum <- jday2[ii] } } cdatxx <- cdatxx[pkthin,] jday2 <- jday2[pkthin] sday2 <- sday2[pkthin] # remove values with missing streamflow anomalies qchkna <- is.na(qanom30[jday2]) | is.na(qanom1[jday2]) cdatxx <- cdatxx[!qchkna,] jday2 <- jday2[!qchkna] sday2 <- sday2[!qchkna] # # create dataframes for input to swaveqexFit # nobs <- length(jday2) icen2xx <- as.character(cdatxx[,2]) icen2 <- rep("u",nobs) icen2[icen2xx=="<"] <- "c" cdatxx[,2] <- as.character(icen2) datyyy <- data.frame(cdatxx[,1],jday2,sday2,log10(cdatxx[,3]),cdatxx[,2],qanom30[jday2],qanom1[jday2]) datyyy[,1] <- as.character(datyyy[,1]) datyyy[,5] <- as.character(datyyy[,5]) names(datyyy) <- c('date','jday','sday',"lgconc","icen",'qanom30','qanom1') cenyr <- nall/365/2 icenall <- rep("u",nall) icenall[jday2[datyyy[,"icen"]=="c"]] <- "c" datsub <- data.frame(datyyy[,1:5],datyyy[,4],rep(1,nobs),rep(NA,nobs),datyyy[,6:7],datyyy[,2]/365-cenyr) datall <- data.frame(datxxx[,1:3],rep(NA,nall),icenall,rep(NA,nall),rep(1,nall),rep(NA,nall),datxxx[,4:5],datxxx[,2]/365-cenyr) datsub[,1] <- as.character(datsub[,1]) datsub[,5] <- as.character(datsub[,5]) datall[,1] <- as.character(datall[,1]) datall[,5] <- as.character(datall[,5]) names(datsub) <- c('date','jday','sday','yobs','icen','yimp','int','swave','qanom30','qanom1','tnd') names(datall) <- c('date','jday','sday','yobs','icen','yimp','int','swave','qanom30','qanom1','tnd') # produce rough plot of observations and flow anomalies xxx <- yrbeg+datsub[,2]/365 yyy <- datsub[,4] icen <- as.character(datsub[,5]) plot(xxx,yyy,xlim=c(yrbeg,yrend+1),type="p",pch=" ") mtext(side=3,line=1,cex=.8,"Observations") for (i in yrbeg:yrend) {abline(v=i) } points(xxx[icen=="u"],yyy[icen=="u"],pch=16,cex=.7) points(xxx[icen=="c"],yyy[icen=="c"],pch=1,cex=.7,col=2) nucen <- sum(icen=="u") ntot <- length(icen) pucen <- 100*round(nucen/ntot,2) mtext(side=3,line=0,paste("numobs=",as.character(ntot)," nuncen=",as.character(nucen),"prucen=",as.character(pucen))) xxx <- yrbeg+datall[,2]/365 yyy <- datall[,9] plot(xxx,yyy,xlim=c(yrbeg,yrend+1),type="l") mtext(side=3,line=1,cex=.8,"Mid-term Flow Anomaly") for (i in yrbeg:yrend) {abline(v=i) } yyy <- datall[,10] plot(xxx,yyy,xlim=c(yrbeg,yrend+1),type="l") mtext(side=3,line=1,cex=.8,"Short-term Flow Anomaly") for (i in yrbeg:yrend) {abline(v=i) } sdqamt <- sd(datall[,"qanom30"],na.rm=T) sdqast <- sd(datall[,"qanom1"],na.rm=T) saveout <- list(paste(runname,qwstnum,sep=""),c(yrbeg,yrend),datsub,datall,qdaily,c(ntot,nucen,pucen),c(sdqamt,sdqast),redconc) } ######################################################################################################### swaveqexFit <- function(datin,outfolder,ncs=100,itrans=1,samcov="full",loc=0,stnd="no",regplot="no") { # # Fits the seawaveqex model using input data prepared by swaveqexMerge # Produces plot file with diagnostic plots and text file with conditional simulations # # datin is a list produced by swaveqexMerge # outfolder is a character variable specifying the folder to put the output files # For example, "outAtrazine\\" puts the files in a folder named outAtrazine in # the default directory. The folder needs to be created beforehand. # ncs is the number of conditional simulations to generate # # Use samcov="full" (the default) if samples are collected throughout the year # Use samcov="partial" if samples are collected during a shorter sampling season # Generally, use samcov="full" if the sampling season is 10 months or more, # samcov="partial" if the sampling season is 6 months or less, and best # judgement if the sampling season is between 6 to 10 months # # itrans specifies the transformation to use for concentration: # itrans=1 (the default) uses a log transformation # itrans=2 uses a square-root transformation # itrans=3 uses a cube-root transformation # itrans <- round(itrans,0) if(itrans<=1) {itrans <- 1 } else if (itrans>3) {itrans <- 3} chtrans <- c("base-10 logarithm","square root","cube root")[itrans] # loc is a nonnegative numeric constant specifying an optional limit of concern # loc=0 (the default) does not include the limit of concern # loc=10 (for example) specifies that the limit of concern is 10 micrograms/per liter if(loc<0) {loc=0} # stnd ("show trend line") specifies whether or not to show the trend line on the # output plots # stnd="no" (the default) does not show the trend line # stnd="yes" shows the trend line # If regplot="yes", a plot is produced showing the fitted values from # the regression model. This plot is comparable to the results from the # seawave-Q model # mxalph is a positive integer that control the maximum allowable value for the parameter alpha # if mxalph=k, where k=1,2, ..., then alpha must be less than the minimum of 1.5 and cswave/(k*sigma) # (Note: mxalph=2 is used for this version) mxalph <- 2 qwstnum <- datin[[1]] yrbeg <- datin[[2]][1] yrend <- datin[[2]][2] datsub <- datin[[3]] datall <- datin[[4]] redconc <- datin[[8]] ### if(redconc>0) { chtrans <- paste(chtrans,"(*concentration reduced by",as.character(redconc),"micrograms per liter)") } else if(redconc<0) { chtrans <- paste(chtrans,"(*concentration increased by",as.character(-redconc),"micrograms per liter)") } ## if itrans>1, over-ride the log transformation if (itrans>1) { datsub[,4] <- (10^datsub[,4])^(1/itrans) datsub[,6] <- (10^datsub[,6])^(1/itrans) datall[,4] <- (10^datall[,4])^(1/itrans) datall[,6] <- (10^datall[,6])^(1/itrans) } qdaily <- datin[[5]] qdate <- as.character(datall[,1]) pdf(file=paste(outfolder,"Plots",qwstnum,".pdf",sep=""),height=6,width=10) par(omi=c(.2,.2,.7,.2)) rlen <- yrend-yrbeg+1 nall <- rlen*365 svcounts <- datin[[6]] svsdflowanoms <- datin[[7]] # # Seasonal wave selection and parameter estimation # svstuff <- swaveqexPESTpdo(datsub,samcov,itrans,mxalph) # Generate conditional simulations # ngen is the number of conditional simulations # parxxx <- svstuff[1:15] spmodcls <- svstuff[16] svpvals <- svstuff[17:20] svstuff <- svstuff[1:15] ngen <- ncs if(ngen>250) {ngen <- 250} tmpout <- swaveqexCSIM(ngen,parxxx,datsub,datall,spmodcls) estannmax <- round(tmpout[[1]],5) p10annmax <- round(tmpout[[5]],5) p90annmax <- round(tmpout[[6]],5) flagnas <- tmpout[[7]] datsub <- data.frame(datsub,tmpout[[2]]) datall <- data.frame(datall,tmpout[[3]]) ##### itrans if(itrans==1) { svcondsim <- round(10^tmpout[[4]],5) } else { estannmax[estannmax<1] <- 1 p10annmax[p10annmax<1] <- 1 p90annmax[p90annmax<1] <- 1 estannmax <- log10(estannmax)^itrans p10annmax <- log10(p10annmax)^itrans p90annmax <- log10(p90annmax)^itrans svcondsim <- tmpout[[4]] svcondsim[svcondsim<0] <- 0 svcondsim <- round(svcondsim^itrans,5) } for (jjj in 1:ngen) { ## add redconc back in svcondsim[,jjj] <- svcondsim[,jjj]+redconc pickxxx <- svcondsim[,jjj] < 0.00001 svcondsim[pickxxx,jjj] <- 0.00001 pickxxx <- is.na(svcondsim[,jjj]) svcondsim[pickxxx,jjj] <- (-0.009) } svcenout <- rep(" ",nall) svcenout[datsub[,"jday"]] <- datsub[,"icen"] svcenout[svcenout=="c"] <- "<" svcenout[svcenout=="u"] <- "o" ##### itrans if(itrans==1) { concorig <- (round(10^datsub[,"yobs"],5)) } else { concorig <- (round(datsub[,"yobs"]^itrans,5)) } ### compute various statistics for observed data jdxx <- datsub[,"jday"] nxx <- yrend-yrbeg+1 nobsxx <- rep(NA,nxx) ssavexx <- rep(NA,nxx); ssminxx <- rep(NA,nxx); ssmaxxx <- rep(NA,nxx) cmaxrecxx <- rep(NA,nxx) for (kxx in 1:nxx) { jdxx1 <- (kxx-1)*365+1 jdxx2 <- kxx*365 pkxx <- jdxx>=jdxx1 & jdxx<=jdxx2 nobsxx[kxx] <- sum(pkxx) if(sum(pkxx)>=2) { jtmpxx <- jdxx[pkxx] ctmpxx <- concorig[pkxx] ntmpxx <- sum(pkxx) jdifxx <- jtmpxx[2:ntmpxx]-jtmpxx[1:(ntmpxx-1)] ssavexx[kxx] <- round(mean(jdifxx),1) ssminxx[kxx] <- min(jdifxx) ssmaxxx[kxx] <- max(jdifxx) cmaxrecxx[kxx] <- max(ctmpxx) } } #### # save output stats in a file # yrxxx <- yrbeg:yrend savestatsxx <- data.frame(yrxxx,nobsxx,ssavexx,ssminxx,ssmaxxx, cmaxrecxx+redconc,estannmax+redconc,p10annmax+redconc,p90annmax+redconc,stringsAsFactors=F) for (jjj in 6:9 ) { savestatsxx[,jjj] <- format(savestatsxx[,jjj],scientific=F) } names(savestatsxx) <- c("year","nobs","ssave","ssmin","ssmax","obsmax","estmax","p10max","p90max") write.table(format(savestatsxx),file=paste(outfolder,"STATS",qwstnum,".txt",sep=""),row.names=F,sep="\t",quote=F) svconcout <- rep(" ",nall) ### svconcout[datsub[,"jday"]] <- as.character(concorig+redconc) qdate <- as.character(datall[,1]) qyear <- as.numeric(substring(qdate,1,4)) qsday <- datall[,3] qout <- round(qdaily,1) qout[is.na(qout)] <- (-9.0) ### itrans if(itrans==1) { muout <- round(10^datall[,"mu"],5) mucondout <- round(10^datall[,"condmean"],5) } else { muout <- datall[,"mu"] ; muout[muout<0] <- 0 muout <- round(muout^itrans,5) mucondout <- datall[,"condmean"] ; mucondout[mucondout<0] <- 0 mucondout <- round(mucondout^itrans,5) } ## ## add back redconc muout <- muout+redconc mucondout <- mucondout+redconc muout[muout<0.00001] <- 0.00001 mucondout[mucondout<0.00001] <- 0.00001 muout[is.na(muout)] <- (-0.009) mucondout[is.na(mucondout)] <- (-0.009) svcondsim <- 1000*svcondsim muout <- 1000*muout mucondout <- 1000*mucondout # # save conditional simulations in a file # savecsim <- data.frame(qdate,qyear,qsday,qout,svconcout,svcenout,muout,mucondout,svcondsim[,1:min(ngen,100)],stringsAsFactors=F) for (jjj in 7:(9+min(ngen,100)-1) ) { savecsim[,jjj] <- format(savecsim[,jjj],scientific=F) } names(savecsim) <- c("date","year","jday","qobs","cobs","crem","estreg","estcmu",paste("csim",as.character(1:min(ngen,100)),sep="")) write.table(format(savecsim),file=paste(outfolder,"CSIMS",qwstnum,".txt",sep=""),row.names=F,sep="\t",quote=F) # # produce plots # hlife <- parxxx[10] ; cmax <- parxxx[8] ; modn <- parxxx[9] swave <- compwaveconvxx(cmax,modn,hlife,spmodcls) sdum <- (1:365)/366 ipk <- floor(360*sdum) ipk[ipk==0] <- 1 swave <- swave[ipk] swaveobs <- swave[datsub[,"sday"]] swaveall <- swave[datall[,"sday"]] datsub[,"swave"] <- swaveobs datall[,"swave"] <- swaveall #### itrans if(itrans==1) { annmax <- log10(estannmax) p10max <- log10(p10annmax) p90max <- log10(p90annmax) } else { annmax <- estannmax^(1/itrans) p10max <- p10annmax^(1/itrans) p90max <- p90annmax^(1/itrans) } condsim <- datall[,"condsim"] condmean <- datall[,"condmean"] yobs <- datsub[,"yobs"] condsimobs <- datsub[,"condsim"] if(itrans>1) { condsim[condsim<0] <- 0 condmean[condmean<0] <- 0 yobs[yobs<0] <- 0 condsimobs[condsimobs<0] <- 0 } yymax <- (max(c(annmax,p90max,condsim,condmean),na.rm=T)) yymin <- floor(min(c(annmax,condsim,condmean),na.rm=T)) if(yymin< (-5)) {yymin <- (-5) } yymax <- ceiling(yymax + (yymax-yymin)*0.25) jdayobs <- datsub[,"jday"] ; dyrobs <- yrbeg+(jdayobs)/365 jdayall <- datall[,"jday"] ; dyrall <- yrbeg+(jdayall)/365 sdayobs <- datsub[,"sday"] ; sdayall <- datsub[,"sday"] muobs <- datsub[,"mu"] muall <- datall[,"mu"] omgobs <- datsub[,"omg"] omgall <- datall[,"omg"] yhatobs <- datsub[,"yhat"]; yhatall <- datall[,"yhat"] yadj2obs <- datsub[,"yadj2"]; yadj2all <- datall[,"yadj2"] yadj3obs <- datsub[,"yadj3"]; yadj3all <- datall[,"yadj3"] yadj4obs <- datsub[,"yadj4"]; yadj4all <- datall[,"yadj4"] regfitvals <- svstuff[1]+svstuff[2]*swaveall + yhatall reguplim <- regfitvals + 2*omgall*svstuff[14] icen <- datsub[,"icen"]=="c" ########## plot1 par(mai=c(.6,1,.1,.2)) plot(dyrall,condsim,type="p",pch=" ",xlim=c(yrbeg-0.2,yrend+1.2),ylim=c(yymin,yymax),xaxt="n",yaxt="n", ylab=" ",xlab=" ",yaxs="i",xaxs="i") for (j in 0:(rlen)) {abline(v=yrbeg+j,col="black",lty=2)} for (j in 0:(rlen-1)) { polygon(c(yrbeg+j,yrbeg+j+1,yrbeg+j+1,yrbeg+j),c(p10max[j+1],p10max[j+1],p90max[j+1],p90max[j+1]), col="light blue",border="medium blue") lines(c(yrbeg+j,yrbeg+j+1),c(annmax[j+1],annmax[j+1]),col="medium blue",lwd=1) if(flagnas[j+1]==T) {text(yrbeg+j+0.5,p90max[j+1]+(yymax-yymin)*0.014,cex=.65,"PR")} mtext(side=1,line=0.5,cex=0.8,at=yrbeg+j+0.5,as.character(yrbeg+j)) } lines(dyrall,condsim,type="l",col="dark grey") points(dyrobs[!icen],condsim[jdayobs[!icen]],pch=16,col="red",cex=0.6) points(dyrobs[icen],condsim[jdayobs[icen]],pch=1,col="red",cex=0.65) if (sum(icen)>0) { abline(h=max(yobs[icen]),col="red") } # add trend line if(stnd=="yes") { xxmin <- yrbeg-0.2; xxmax <- yrend+1.2; xxcen <- (xxmin+xxmax)/2 lines(c(xxmin,xxmax),svstuff[1]+svstuff[5]*c(xxmin-xxcen,xxmax-xxcen),col="medium blue",lty=4,lwd=1.5) } # add line for loc if (loc>0) { loctmp <- loc if(itrans==1) {loctmp <- log10(loc) } else { loctmp <- loc^(1/itrans)} if(loctmp< yymax-0.01*(yymax-yymin)) {abline(h=loctmp,col="red",lty=4,lwd=1.5)} } #### itrans if(yymax-yymin>3) { yytck <- yymin:yymax } else { yytck <- seq(yymin,yymax,0.5) } if(itrans==1) {yylab <- format(10^yytck,big.mark=",",scientific=F,drop0trailing=T) } else { yylab <- format(yytck^itrans,big.mark=",",scientific=F,drop0trailing=T) } for (j in yytck) {lines(c(yrbeg-0.2,yrbeg-0.2+0.01*rlen),c(j,j)) ; lines(c(yrend+1.2-0.01*rlen,yrend+1.2),c(j,j))} mtext(side=2,line=0.5,adj=1,cex=0.8,at=yytck,yylab,las=1) yydel <- (yymax-yymin)*0.025 if (stnd=="yes" & loc>0) {xxl <- 7 } else {xxl <- 6} polygon(c(yrbeg+.16*rlen,yrbeg+.9*rlen,yrbeg+.9*rlen,yrbeg+.16*rlen), c(yymax-xxl*yydel,yymax-xxl*yydel,yymax-yydel,yymax-yydel),col="white") text(yrbeg+.5*rlen,yymax-1.9*yydel,cex=.65,"EXPLANATION") lines(c(yrbeg+.5*rlen,yrbeg+.54*rlen),rep(yymax-3*yydel,2),col="dark grey",lwd=1.5) text(yrbeg+.55*rlen,yymax-3*yydel,cex=.65,adj=0,paste("Conditional trace (first of",as.character(ncs),"traces)")) lines(c(yrbeg+.5*rlen,yrbeg+.54*rlen),rep(yymax-4*yydel,2),col="dark blue",lwd=1.5) text(yrbeg+.55*rlen,yymax-4*yydel,cex=.65,adj=0,"Estimated annual maximum (with 80% error bounds)") if (stnd=="yes") { lines(c(yrbeg+.5*rlen,yrbeg+.54*rlen),rep(yymax-5*yydel,2),col="medium blue",lwd=1.5,lty=4) text(yrbeg+.55*rlen,yymax-5*yydel,cex=.65,adj=0,"Estimated trend") xxl <- 6 } else {xxl <- 5} if (loc>0) { lines(c(yrbeg+.5*rlen,yrbeg+.54*rlen),rep(yymax-xxl*yydel,2),col="red",lwd=1.5,lty=4) text(yrbeg+.55*rlen,yymax-xxl*yydel,cex=.65,adj=0,"Limit of concern") } points(yrbeg+0.21*rlen,yymax-3*yydel,pch=16,col="red",cex=0.6) text(yrbeg+.24*rlen,yymax-3*yydel,cex=.65,adj=0,"Uncensored observation") points(yrbeg+0.21*rlen,yymax-4*yydel,pch=1,col="red",cex=0.65) text(yrbeg+.24*rlen,yymax-4*yydel,cex=.65,adj=0,"Simulated censored observation") lines(c(yrbeg+.19*rlen,yrbeg+.23*rlen),rep(yymax-5*yydel,2),col="red") text(yrbeg+.24*rlen,yymax-5*yydel,cex=.65,adj=0,"Highest censoring limit") mtext(side=1,line=2,cex=0.9,"Calendar year") mtext(side=2,line=3,cex=0.9,"Concentration, in micrograms per liter") mtext(side=3,outer=T,line=2,cex=1.0,qwstnum) mtext(side=3,outer=T,line=1,cex=1.0,paste("Transformation:",chtrans)) if(sum(flagnas)>0) {mtext(side=3,outer=T,line=0,cex=1.0,"PR indicates annual maximum is based on partial record")} ###### plot 1A (seawave-Q regression plot) if (regplot=="yes") { plot(dyrall,regfitvals,type="p",pch=" ",xlim=c(yrbeg-0.2,yrend+1.2),ylim=c(yymin,yymax),xaxt="n",yaxt="n", ylab=" ",xlab=" ",yaxs="i",xaxs="i") for (j in 0:(rlen)) {abline(v=yrbeg+j,col="black",lty=2)} for (j in 0:(rlen-1)) { mtext(side=1,line=0.5,cex=0.8,at=yrbeg+j+0.5,as.character(yrbeg+j)) } points(dyrobs[!icen],yobs[!icen],pch=16,col="red",cex=0.6) points(dyrobs[icen],yobs[icen],pch=1,col="red",cex=0.65) lines(dyrall,regfitvals,type="l",col="black") lines(dyrall,reguplim,type="l",col="dark grey") if (sum(icen)>0) { abline(h=max(yobs[icen]),col="red") } # add trend line if(stnd=="yes") { xxmin <- yrbeg-0.2; xxmax <- yrend+1.2; xxcen <- (xxmin+xxmax)/2 lines(c(xxmin,xxmax),svstuff[1]+svstuff[5]*c(xxmin-xxcen,xxmax-xxcen),col="medium blue",lty=4,lwd=1.5) } # add line for loc if (loc>0) { loctmp <- loc if(itrans==1) {loctmp <- log10(loc) } else { loctmp <- loc^(1/itrans)} if(loctmp< yymax-0.01*(yymax-yymin)) {abline(h=loctmp,col="red",lty=4,lwd=1.5)} } #### itrans if(yymax-yymin>3) { yytck <- yymin:yymax } else { yytck <- seq(yymin,yymax,0.5) } if(itrans==1) {yylab <- format(10^yytck,big.mark=",",scientific=F,drop0trailing=T) } else { yylab <- format(yytck^itrans,big.mark=",",scientific=F,drop0trailing=T) } for (j in yytck) {lines(c(yrbeg-0.2,yrbeg-0.2+0.01*rlen),c(j,j)) ; lines(c(yrend+1.2-0.01*rlen,yrend+1.2),c(j,j))} mtext(side=2,line=0.5,adj=1,cex=0.8,at=yytck,yylab,las=1) yydel <- (yymax-yymin)*0.025 if (stnd=="yes" & loc>0) {xxl <- 7 } else {xxl <- 6} polygon(c(yrbeg+.16*rlen,yrbeg+.9*rlen,yrbeg+.9*rlen,yrbeg+.16*rlen), c(yymax-xxl*yydel,yymax-xxl*yydel,yymax-yydel,yymax-yydel),col="white") text(yrbeg+.5*rlen,yymax-1.9*yydel,cex=.65,"EXPLANATION") lines(c(yrbeg+.5*rlen,yrbeg+.54*rlen),rep(yymax-3*yydel,2),col="black",lwd=1.5) text(yrbeg+.55*rlen,yymax-3*yydel,cex=.65,adj=0,"Fitted values from regression model") lines(c(yrbeg+.5*rlen,yrbeg+.54*rlen),rep(yymax-4*yydel,2),col="dark grey",lwd=1.5) text(yrbeg+.55*rlen,yymax-4*yydel,cex=.65,adj=0,"Upper 95% prediction limits from regression model") if (stnd=="yes") { lines(c(yrbeg+.5*rlen,yrbeg+.54*rlen),rep(yymax-5*yydel,2),col="medium blue",lwd=1.5,lty=4) text(yrbeg+.55*rlen,yymax-5*yydel,cex=.65,adj=0,"Estimated trend") xxl <- 6 } else {xxl <- 5} if (loc>0) { lines(c(yrbeg+.5*rlen,yrbeg+.54*rlen),rep(yymax-xxl*yydel,2),col="red",lwd=1.5,lty=4) text(yrbeg+.55*rlen,yymax-xxl*yydel,cex=.65,adj=0,"Limit of concern") } points(yrbeg+0.21*rlen,yymax-3*yydel,pch=16,col="red",cex=0.6) text(yrbeg+.24*rlen,yymax-3*yydel,cex=.65,adj=0,"Uncensored observation") points(yrbeg+0.21*rlen,yymax-4*yydel,pch=1,col="red",cex=0.65) text(yrbeg+.24*rlen,yymax-4*yydel,cex=.65,adj=0,"Censored observation") lines(c(yrbeg+.19*rlen,yrbeg+.23*rlen),rep(yymax-5*yydel,2),col="red") text(yrbeg+.24*rlen,yymax-5*yydel,cex=.65,adj=0,"Highest censoring limit") mtext(side=1,line=2,cex=0.9,"Calendar year") mtext(side=2,line=3,cex=0.9,"Concentration, in micrograms per liter") mtext(side=3,outer=T,line=2,cex=1.0,qwstnum) mtext(side=3,outer=T,line=1,cex=1.0,paste("Transformation:",chtrans)) if(sum(flagnas)>0) {mtext(side=3,outer=T,line=0,cex=1.0,"PR indicates annual maximum is based on partial record")} } ########### plot2 par(mai=c(.6,1.5,.5,.7)) yy <- condsimobs-yhatobs plot(sdayobs,yy,type="p",pch=" ",xlim=c(-15,380),ylim=c(yymin,yymax),xaxt="n",yaxt="n",ylab=" ",xlab=" ",yaxs="i",xaxs="i") for (j in 1:365) { if(swaveall[j]>0) {polygon(c(j,j+1,j+1,j),c(yymin+.2*yydel,yymin+.2*yydel,yymax-.2*yydel,yymax-.2*yydel),col="light gray",border=NA)} } points(sdayobs[!icen],yy[!icen],pch=16,col="red",cex=0.6) points(sdayobs[icen],yy[icen],pch=1,col="red",cex=0.65) lines(1:365,svstuff[1]+svstuff[2]*swaveall[1:365],lwd=1.5) lines(1:365,svstuff[1]+svstuff[2]*swaveall[1:365]+2*omgall[1:365]*svstuff[14],lwd=1.5,lty=4) lines(1:365,svstuff[1]+svstuff[2]*swaveall[1:365]-2*omgall[1:365]*svstuff[14],lwd=1.5,lty=4) for (j in yytck) {lines(c(-15,-10),c(j,j)) ; lines(c(375,380),c(j,j))} mtext(side=2,line=0.5,adj=1,cex=0.8,at=yytck,yylab,las=1) mocum <- c(0,31,59,90,120,151,181,212,243,273,304,334) for (j in c(mocum,365)) { abline(v=j,col="grey",lty=2) } mtext(side=1,line=0.5,cex=0.8,at=mocum+15,c("Jan","Feb","Mar","Apr","May","June","July","Aug","Sept","Oct","Nov","Dec")) polygon(c(.61*365,370,370,.61*365), c(yymax-8*yydel,yymax-8*yydel,yymax-yydel,yymax-yydel),col="white") text(.75*365,yymax-2*yydel,cex=.6,"EXPLANATION") lines(c(.62*365,.68*365),rep(yymax-3*yydel,2),lty=1,lwd=1.5) text(.7*365,yymax-3*yydel,cex=.6,adj=0,"Seasonal wave") lines(c(.62*365,.68*365),rep(yymax-4*yydel,2),lty=4,lwd=1.5) text(.7*365,yymax-4*yydel,cex=.6,adj=0,"Seasonal wave plus and minus two times SSD") points(0.65*365,yymax-5*yydel,pch=16,col="red",cex=0.6) text(.7*365,yymax-5*yydel,cex=.6,adj=0,"Uncensored observation") points(0.65*365,yymax-6*yydel,pch=1,col="red",cex=0.65) text(.7*365,yymax-6*yydel,cex=.6,adj=0,"Simulated censored observation") text(.7*365,yymax-7*yydel,cex=.6,adj=0,"SSD, seasonal standard deviation") # text(.7*365,yymax-8*yydel,cex=.6,adj=0,"Shaded region, high conc. season(s)") mtext(side=2,line=3,cex=0.9,"Adjusted concentration, in micrograms per liter") mtext(side=3,outer=T,line=2,cex=1.0,qwstnum) mtext(side=3,outer=T,line=1,cex=1.0,paste("Transformation:",chtrans)) mtext(side=3,outer=T,line=0,cex=1,"Adjusted concentration (minus trend and flow-related variability) versus month") mtext(side=3,outer=T,line=-1,cex=1,"Shaded region indicates high concentration season(s)") ###### plot3 par(mai=c(.6,1.5,.5,.7)) yy <- condsimobs-yadj2obs xx2 <- datsub[,"qanom30"] xx2min <- floor(min(xx2)*2)/2; xx2max <- ceiling(max(xx2)*2)/2 xx2del <- (xx2max-xx2min)*0.01 ; yydel <- (yymax-yymin)*0.015 plot(xx2,yy,type="p",pch=" ",xlim=c(xx2min,xx2max),ylim=c(yymin,yymax),xaxt="n",yaxt="n",ylab=" ",xlab=" ",yaxs="i",xaxs="i") points(xx2[!icen],yy[!icen],pch=16,col="red",cex=0.6) points(xx2[icen],yy[icen],pch=1,col="red",cex=0.65) abline(svstuff[1],svstuff[3]) for (j in yytck) {lines(c(xx2min,xx2min+xx2del),c(j,j)) ; lines(c(xx2max-xx2del,xx2max),c(j,j))} mtext(side=2,line=0.5,adj=1,cex=0.8,at=yytck,yylab,las=1) for (j in seq(xx2min,xx2max,0.5)) { mtext(side=1,line=0.5,cex=0.8,at=j,as.character(j)) lines(c(j,j),c(yymin,yymin+yydel)) ; lines(c(j,j),c(yymax-yydel,yymax)) } mtext(side=2,line=3,cex=0.9,"Adjusted concentration, in micrograms per liter") mtext(side=1,line=2,cex=0.9,"Mid-term flow anomaly") mtext(side=3,outer=T,line=2,cex=1.0,qwstnum) mtext(side=3,outer=T,line=1,cex=1.0,paste("Transformation:",chtrans)) mtext(side=3,outer=T,line=0,cex=1,"Adjusted concentration (minus seasonality, trend, and STFA) versus MTFA") mtext(side=3,outer=T,line=-1,cex=1,paste("(coef=", as.character(round(svstuff[3],3)),", pval=",as.character(round(svpvals[2],3)),")")) ###### plot4 yy <- condsimobs-yadj3obs xx2 <- datsub[,"qanom1"] xx2min <- floor(min(xx2)*2)/2; xx2max <- ceiling(max(xx2)*2)/2 xx2del <- (xx2max-xx2min)*0.01 ; yydel <- (yymax-yymin)*0.015 plot(xx2,yy,type="p",pch=" ",xlim=c(xx2min,xx2max),ylim=c(yymin,yymax),xaxt="n",yaxt="n",ylab=" ",xlab=" ",yaxs="i",xaxs="i") points(xx2[!icen],yy[!icen],pch=16,col="red",cex=0.6) points(xx2[icen],yy[icen],pch=1,col="red",cex=0.65) abline(svstuff[1],svstuff[4]) for (j in yytck) {lines(c(xx2min,xx2min+xx2del),c(j,j)) ; lines(c(xx2max-xx2del,xx2max),c(j,j))} mtext(side=2,line=0.5,adj=1,cex=0.8,at=yytck,yylab,las=1) for (j in seq(xx2min,xx2max,0.5)) { mtext(side=1,line=0.5,cex=0.8,at=j,as.character(j)) lines(c(j,j),c(yymin,yymin+yydel)) ; lines(c(j,j),c(yymax-yydel,yymax)) } mtext(side=2,line=3,cex=0.9,"Adjusted concentration, in micrograms per liter") mtext(side=1,line=2,cex=0.9,"Short-term flow anomaly") mtext(side=3,outer=T,line=2,cex=1.0,qwstnum) mtext(side=3,outer=T,line=1,cex=1.0,paste("Transformation:",chtrans)) mtext(side=3,outer=T,line=0,cex=1,"Adjusted concentration (minus seasonality, trend, and MTFA) versus STFA") mtext(side=3,outer=T,line=-1,cex=1,paste("(coef=", as.character(round(svstuff[4],3)),", pval=",as.character(round(svpvals[3],3)),")")) ###### plot4a yy <- condsimobs-yadj4obs xx2 <- dyrobs xx2min <- floor(min(xx2)); xx2max <- ceiling(max(xx2)); xx2cen <- (xx2max+xx2min)/2 xx2del <- (xx2max-xx2min)*0.01 ; yydel <- (yymax-yymin)*0.015 plot(xx2,yy,type="p",pch=" ",xlim=c(xx2min,xx2max),ylim=c(yymin,yymax),xaxt="n",yaxt="n",ylab=" ",xlab=" ",yaxs="i",xaxs="i") points(xx2[!icen],yy[!icen],pch=16,col="red",cex=0.6) points(xx2[icen],yy[icen],pch=1,col="red",cex=0.65) lines(c(xx2min,xx2max),svstuff[1]+svstuff[5]*c(xx2min-xx2cen,xx2max-xx2cen)) for (j in yytck) {lines(c(xx2min,xx2min+xx2del),c(j,j)) ; lines(c(xx2max-xx2del,xx2max),c(j,j))} mtext(side=2,line=0.5,adj=1,cex=0.8,at=yytck,yylab,las=1) for (j in 0:(rlen-1)) { if(j<(rlen-1)) {abline(v=yrbeg+j+1,col="black",lty=2)} mtext(side=1,line=0.5,cex=0.8,at=yrbeg+j+0.5,as.character(yrbeg+j)) } mtext(side=2,line=3,cex=0.9,"Adjusted concentration, in micrograms per liter") mtext(side=1,line=2,cex=0.9,"Calendar year") mtext(side=3,outer=T,line=2,cex=1.0,qwstnum) mtext(side=3,outer=T,line=1,cex=1.0,paste("Transformation:",chtrans)) mtext(side=3,outer=T,line=0,cex=1,"Adjusted concentration (minus seasonality and flow-related variablility) versus year") mtext(side=3,outer=T,line=-1,cex=1,paste("(coef=", as.character(round(svstuff[5],3)),", pval=",as.character(round(svpvals[4],3)),")")) ###### plot5 yy <- (condsimobs-muobs)/omgobs/svstuff[14] yy[yy>3] <- 3 yy[yy< (-3)] <- (-3) plot(sdayobs,yy,type="p",pch=" ",xlim=c(-15,380),ylim=c(-3.2,3.2),xaxt="n",yaxt="n",ylab=" ",xlab=" ",yaxs="i",xaxs="i") for (j in 1:365) { if(swaveall[j]>0) {polygon(c(j,j+1,j+1,j),c(-3,-3,3,3),col="light gray",border=NA)} } points(sdayobs[!icen],yy[!icen],pch=16,col="red",cex=0.6) points(sdayobs[icen],yy[icen],pch=1,col="red",cex=0.65) for (j in (-3):3) { mtext(side=2,line=0.5,adj=1,cex=0.8,at=j,as.character(j),las=1) if (abs(j)<=3) {abline(h=j,col="grey",lty=2) } } abline(h=0) for (j in c(mocum,365)) { abline(v=j,col="grey",lty=2) } mtext(side=1,line=0.5,cex=0.8,at=mocum+15,c("Jan","Feb","Mar","Apr","May","June","July","Aug","Sept","Oct","Nov","Dec")) mtext(side=2,line=3,cex=0.9,"Normalized residual (dimensionless)") mtext(side=3,outer=T,line=2,cex=1.0,qwstnum) mtext(side=3,outer=T,line=1,cex=1.0,paste("Transformation:",chtrans)) mtext(side=3,outer=T,line=0,cex=1,"Normalized residuals versus month") mtext(side=3,outer=T,line=-1,cex=1,"Shaded region indicates high concentration season(s)") plot(dyrobs,yy,type="p",pch=" ",xlim=c(yrbeg-0.2,yrend+1.2),ylim=c(-3.2,3.2),xaxt="n",yaxt="n",ylab=" ",xlab=" ",yaxs="i",xaxs="i") points(dyrobs[!icen],yy[!icen],pch=16,col="red",cex=0.6) points(dyrobs[icen],yy[icen],pch=1,col="red",cex=0.7) for (j in (-3):3) { mtext(side=2,line=0.5,adj=1,cex=0.8,at=j,as.character(j),las=1) if (abs(j)<=3) {abline(h=j,col="grey",lty=2) } } for (j in 0:(rlen-1)) { abline(v=yrbeg+j+1,col="grey",lty=2) mtext(side=1,line=0.5,cex=0.8,at=yrbeg+j+0.5,as.character(yrbeg+j)) abline(v=yrbeg,col="grey",lty=2) } mtext(side=2,line=3,cex=0.9,"Normalized residual (dimensionless)") mtext(side=3,outer=T,line=2,cex=1.0,qwstnum) mtext(side=3,outer=T,line=1,cex=1.0,paste("Transformation:",chtrans)) mtext(side=3,outer=T,line=0,cex=1,"Normalized residuals versus year") ###### correlogram plot rxx <- yy/sd(yy) txx <- dyrobs nxx <- length(rxx) nxx2 <- nxx*(nxx-1)/2 dtmp <- rep(NA,nxx2) vtmp <- rep(NA,nxx2) cxx <- 0 for(ixx in 1:(nxx-1)) { for(jxx in (ixx+1):nxx) { cxx <- cxx+1 dtmp[cxx] <- txx[jxx]-txx[ixx] vtmp[cxx] <- (rxx[jxx]-rxx[ixx])^2 } } dtmp <- dtmp*365 pcktmp <- dtmp<60 vtmp <- vtmp[pcktmp] dtmp <- dtmp[pcktmp] evxx <- seq(5,55,5) ; nxx <- length(evxx) ; hwid <- rep(5,nxx) evyy <- rep(NA,nxx) ; evyy10 <- evyy; evyy90 <- evyy for (j in 1:nxx) { evpck <- dtmp>=(evxx[j]-hwid[j]) & dtmp<=(evxx[j]+hwid[j]) if(sum(evpck)>=10) { nsub <- sum(evpck); vsub <- vtmp[evpck]^1 ; vrep <- rep(NA,1000) for (irep in 1:1000) { pckrep <- floor(runif(nsub)*nsub)+1 vrep[irep] <- mean(vsub[pckrep]) } evyy10[j] <- quantile(vrep,p=0.025) evyy90[j] <- quantile(vrep,p=0.975) evyy[j] <- mean(vsub) } } evyy <- 1 - 0.5*evyy ; evyy10 <- 1-0.5*evyy10; evyy90 <- 1-0.5*evyy90 rhoxx <- parxxx[13] ylow <- -0.4 ; yup <- 1.2 ytck <- .02*(yup-ylow) xlow <- 0; xup <- 60; xtck <- .012*(xup-xlow) par(mai=c(0.8,2,.3,1.2)) plot(evxx,evyy,pch=' ',cex=.8,xlim=c(xlow,xup),ylim=c(ylow,yup),xaxs='i',yaxs='i',xaxt='n',yaxt='n',xlab='',ylab='') for (j in seq(xlow,xup,10)) { mtext(side=1,line=.5,at=j,cex=.8,as.character(j)) lines(c(j,j),c(ylow,ylow+ytck),lwd=1); lines(c(j,j),c(yup,yup-ytck),lwd=1) } for (j in seq(ylow,yup,.2)) {mtext(side=2,line=.5,at=j,adj=1,las=1,cex=.8,as.character(j)) lines(y=c(j,j),x=c(xlow,xlow+xtck),lwd=1); lines(y=c(j,j),x=c(xup,xup-xtck),lwd=1) } for (j in 1:nxx) { polygon(c(evxx[j]-0.5,evxx[j]+0.5,evxx[j]+0.5,evxx[j]-0.5),c(evyy10[j],evyy10[j],evyy90[j],evyy90[j]), col="light gray",border="gray") } lines(evxx,evyy,type='p',pch=16,cex=.8,col=1) lines(seq(0,xup,.1),rhoxx^seq(0,xup,.1),lwd=2,col=1) abline(h=0) mtext(side=2,line=3,cex=.9,'Correlation between normalized residuals') mtext(side=1,line=2,cex=.9,'Time between observations, in days') if (mean(ssavexx,na.rm=T)<60) {abline(v=mean(ssavexx,na.rm=T),col=2,lwd=1.5)} polygon(c(31.2,58.6,58.6,31.2),c(.86,.86,1.13,1.13),col="white") text(50-8,1.1,cex=.65,"EXPLANATION") lines(c(40,43)-8,c(1.05,1.05),lwd=1.5) text(44-8,1.05,adj=0,cex=.65,"Fitted exponential correlation function") polygon(c(41,42,42,41)-8,c(.93,.93,1.03,1.03),col="light gray",border="gray") points(41.5-8,0.98,cex=.8,pch=16) text(44-8,.98,adj=0,cex=.65,"Empirical correlogram (with 95-pct conf. limits)") lines(c(40,43)-8,c(.9,.9),col=2,lwd=1.5) text(44-8,0.9,adj=0,cex=0.65,lwd=1.5,"Average time between successive observations") mtext(side=3,outer=T,line=2,cex=1.0,qwstnum) mtext(side=3,outer=T,line=1,cex=1.0,paste("Transformation:",chtrans)) chcts <- as.character(round((-1/log(rhoxx)),1)) mtext(side=3,outer=T,line=0,cex=1,"Correlation function of normalized residuals") mtext(side=3,outer=T,line=-1,cex=1,paste("(estimated CTS = ",chcts," days)",sep="")) dev.off() svstuff[13] <- -1/log(svstuff[13]) svstuff <- round(svstuff,4)[c(1:5,8:10,12:15)] svstuff <- c(svstuff[1:5],spmodcls,svstuff[c(7,8,6,11,9,10,12)]) svstuff <- c(svstuff,svcounts,svpvals,svsdflowanoms) svstuff <- c(svstuff[1:20],yrbeg,yrend,rlen,svstuff[21:22]) svnames <- c("int","cswave","cmtfa","cstfa","ctnd","wmcls","wmodno","hlife","wshft","sigma","alph","cts","n2LLIK", "nobs","nucen","prucen","pswave","pmtfa","pstfa","ptnd","yrbeg","yrend","rlen","sdmtfa","sdstfa") svstuff <- svstuff[c(21,22,23,14,15,16,1,2,17,3,18,4,19,5,20,6:13,24,25)] svnames <- svnames[c(21,22,23,14,15,16,1,2,17,3,18,4,19,5,20,6:13,24,25)] list(qwstnum,svnames,svstuff) } ######################################################################################################### swaveqexPESTpdo <- function(dat,samcov,itrans,mxalph) { # # Parameter estimation for seawaveqex model # clog <- dat[,"yobs"] centmp <- dat[,"icen"]=="c" # # The seasonal wave is computed using a period of 360 "days" # rather than 365 because it is required to have the period # be a multiple of 12. The following 3 lines define variables # that are used to "stretch" the period out to 365 days # tseas <- dat[,"sday"]/366 ipk <- floor(360*tseas) ipk[ipk==0] <- 1 # # Populate variables for julian day, mid-term and short-term # flow anomalies. and linear trend # jd <- dat[,"jday"] qamt <- dat[,"qanom30"] qast <- dat[,"qanom1"] lintnd <- dat[,"tnd"] # # Selection of the best seasonal wave model # The loop indeces correspond to the following cases: # - For wave class 1 (spmodcls=1), j is the choice of pulse input model # and k is the choice of the "half-life" (30, 60, 90, or 120 days) # - For wave class 2 (spmodcls=2) and k=1 or 2, j is the choice of pulse # input model from subclass A and k is the choice of "half-life" (30 or 60 days) # - For wave class 2 and k=3 or 4, j is the choice of pulse input model # from subclass B and k-2 is the choice of the "half-life" (30 or 60 days) #### find best wave from class 1 spmodcls <- 1 if(samcov=="partial") {kmax <- 2 }else {kmax <- 4} parxxx <- matrix(nrow=24,ncol=7) pvalxxx <- matrix(nrow=24,ncol=4) svcmaxt <- rep(NA,24) for (j in 1:6) { for (k in 1:kmax) { j2 <- (j-1)*4+k indcen <- !centmp # # For the specified wave model, the phase shift is selected from # 24 choices to maximize the likelihood function liksave <- rep(NA,24) for (imx in 1:24) { cmaxt <- imx/24 awave <- compwaveconvxx(cmaxt,j,k,spmodcls) wavest <- awave[ipk] tmpouta <- survreg(Surv(time=clog,time2=indcen,type="left") ~ wavest+qamt+qast+lintnd,dist="gaussian") liksave[imx] <- -tmpouta$loglik[2] #### restrict to positive wave coefficients if (itrans==1) { minint <- (-3) } else if (itrans==2) { minint <- .001^0.5 } else if (itrans==3) { minint <- .001^(1/3) } chksign <- summary(tmpouta)$table[2,1] chkint <- summary(tmpouta)$table[1,1] if(chksign<0 | chkint< minint) {liksave[imx] <- NA} } pckone <- order(liksave)[1] svcmaxt[j2] <- pckone/24 cmaxt <- pckone/24 awave <- compwaveconvxx(cmaxt,j,k,spmodcls) wavest <- awave[ipk] # # Now save the output for the wave model with the best phase shift # tmpouta <- survreg(Surv(time=clog,time2=indcen,type="left") ~ wavest+qamt+qast+lintnd,dist="gaussian") parxxx[j2,] <- c(summary(tmpouta)$table[1:5,1],tmpouta$scale,-tmpouta$loglik[2]) pvalxxx[j2,] <- c(summary(tmpouta)$table[2:5,4]) } } #### find best wave from class 2 spmodcls <- 2 parxxx2 <- matrix(nrow=24,ncol=7) pvalxxx2 <- matrix(nrow=24,ncol=4) svcmaxt2 <- rep(NA,24) for (j in 1:6) { for (k in 1:4) { j2 <- (j-1)*4+k indcen <- !centmp # # For the specified wave model, the phase shift is selected from # 24 choices to maximize the likelihood function liksave2 <- rep(NA,24) for (imx in 1:24) { cmaxt <- imx/24 awave <- compwaveconvxx(cmaxt,j,k,spmodcls) wavest <- awave[ipk] tmpouta <- survreg(Surv(time=clog,time2=indcen,type="left") ~ wavest+qamt+qast+lintnd,dist="gaussian") liksave2[imx] <- -tmpouta$loglik[2] #### restrict to positive wave coefficients chksign <- summary(tmpouta)$table[2,1] chkint <- summary(tmpouta)$table[1,1] if(chksign<0 | chkint< minint) {liksave[imx] <- NA} } pckone2 <- order(liksave2)[1] svcmaxt2[j2] <- pckone2/24 cmaxt2 <- pckone2/24 awave <- compwaveconvxx(cmaxt2,j,k,spmodcls) wavest <- awave[ipk] # # Now save the output for the wave model with the best phase shift # tmpouta <- survreg(Surv(time=clog,time2=indcen,type="left") ~ wavest+qamt+qast+lintnd,dist="gaussian") parxxx2[j2,] <- c(summary(tmpouta)$table[1:5,1],tmpouta$scale,-tmpouta$loglik[2]) pvalxxx2[j2,] <- c(summary(tmpouta)$table[2:5,4]) } } # Select the best wave model for each class. If the regression coefficient for # the seasonal wave is negative do not choose that model because # the coefficient should be positive to make sense. Penalize class 2 liklihood # by adding 2. Compare likelihoods and choose best class. # likall <- (parxxx[,7]) ; likall[parxxx[,2]<0] <- NA pckone <- order(likall)[1] likval1 <- likall[pckone] likall2 <- (parxxx2[,7]) ; likall2[parxxx2[,2]<0] <- NA pckone2 <- order(likall2)[1] likval2 <- likall2[pckone2]+0 if ((likval1<=likval2 & samcov=="full") | samcov=="partial") { mod1 <- floor((pckone-1)/4)+1 hlife1 <- pckone-(mod1-1)*4 cmaxt <- svcmaxt[pckone] svpars <- c(parxxx[pckone,],cmaxt,mod1,hlife1) svpvals <- round(pvalxxx[pckone,],3) spmodcls <- 1 } else if ((likval1 > likval2 & samcov=="full")){ mod1 <- floor((pckone2-1)/4)+1 hlife1 <- pckone2-(mod1-1)*4 cmaxt <- svcmaxt2[pckone2] svpars <- c(parxxx2[pckone2,],cmaxt,mod1,hlife1) svpvals <- round(pvalxxx2[pckone2,],3) spmodcls <- 2 } awave <- compwaveconvxx(cmaxt,mod1,hlife1,spmodcls) wavest <- awave[ipk] wtmp <- .5+wavest # # Now we have the best wave model and the estimated regression coefficients. # We are ready to compute the regression residuals (ytil) and estimate sigma, # alpha, and delta. Note that rho is the lag-1 correlation coefficient and # delta (the correlation time scale) is -1/ln(rho) # ytil <- (clog-svpars[1]-svpars[2]*wavest-svpars[3]*qamt-svpars[4]*qast-svpars[5]*lintnd) icen <- !indcen # # Set starting values for rho sigma and alpha # Disregard the variable nrat. It is left-over from a previous version and is # always equal to zero nrat <- 0.0 rho0 <- 0.9 alph0 <- 0 sig0 <- svpars[6]/(1+nrat)^.5 wcxxx0 <- svpars[2] maxalphxx <- min(1.5,abs(wcxxx0)/(mxalph*sig0)) # # Now find maximum (pseudo) likelihood estimates of sigma, alpha, and rho # tmpout <- nlminb(objective=evalmodlikxx,control=list(rel.tol=0.0001,x.tol=0.001,iter.max=50),start=c(alph0,log(rho0)), lower=c(0,log(.7515)),upper=c(maxalphxx,log(.9835)),sig0=sig0, nrat=nrat, ytil=ytil, jd=jd, icen=icen, wavest=wavest, ifixs=0) rho1 <- exp(tmpout$par[2]) alph1 <- tmpout$par[1] sig1 <- estsigxx(rho1,nrat,sig0,ytil,jd,icen,wavest,alph1) alph1 <- min(alph1,abs(wcxxx0)/(mxalph*sig1)) pars1 <- c(alph1,log(rho1)) stlik1 <- evalmodlikxx(pars1,nrat,sig1,ytil,jd,icen,wavest,1) # # Ready to save the parameter estimates etc # svout <- round(c(svpars,alph0,alph1,rho1,sig1,stlik1,spmodcls,svpvals),4) svout } ######################################################################################################### estsigxx <- function(rho,nrat,sig0,ytil,jd,icen,wavest,alph) { omg <- (1+alph*(0.0+wavest)) del <- (1+nrat/omg^2)^0.5 y <- ytil/omg/del ntot <- length(y) ncen <- sum(icen) sig1 <- sig0 for (jit in 1:20) { if(!is.na(sig1) & (jit==1 | abs(sig1-sig0)/sig0 > .005)) { sig0 <- sig1 neff <- 0 ssres <- 0 yimp <- y if(icen[1]) { yimp[1] <- sig0*mtmvnorm(mean=0,sigma=1,upper=y[1]/sig0)[[1]] } for (i in 2:ntot) { cxx <- rho^(jd[i]-jd[i-1])/del[i]/del[i-1] if (!icen[i]) { ztmp <- (y[i]-cxx*yimp[i-1])/(1-cxx^2)^0.5/sig0 if(ztmp>2.5) {ztmp <- 2.5 } else if (ztmp < (-2.5)) {ztmp <- (-2.5) } ssres <- ssres+ztmp^2*sig0^2 neff <- neff+1 yimp[i] <- ztmp*(1-cxx^2)^0.5*sig0 + cxx*yimp[i-1] } else { ztmp <- (y[i]-cxx*yimp[i-1])/(1-cxx^2)^0.5/sig0 if(ztmp<(-2.5)) {ztmp <- (-2.5) } neff <- neff-sig0^2*dnorm(ztmp)/pnorm(ztmp)*ztmp yimp[i] <- mtmvnorm(mean=0,sigma=1,upper=ztmp)[[1]] yimp[i] <- yimp[i]*(1-cxx^2)^0.5*sig0 + cxx*yimp[i-1] } } sig1 <- (ssres/neff)^0.5 # if (is.na(sig1)) {sig1 <- sig0} ###### if (is.na(sig1)) {sig1 <- sig0 } else if (sig1>1.05*sig0) {sig1 <- 1.05*sig0} } } sig1 } ######################################################################################################### evalmodlikxx <- function(parin,nrat,sig0,ytil,jd,icen,wavest,ifixs) { #### if ifixs=0, find updated process variance (signew) given rho, nrat, and starting value (sig0) rho <- exp(parin[2]) alph <- parin[1] if (ifixs==0) { signew <- estsigxx(rho,nrat,sig0,ytil,jd,icen,wavest,alph) } else { signew <- sig0 } #### compute conditional mimimum of -2logL omg <- (1+alph*(wavest)) del <- (1+nrat/omg^2)^0.5 y <- ytil/omg/del ntot <- length(y) sig0 <- signew cumlik <- 0 yimp <- y if(icen[1]) { yimp[1] <- sig0*mtmvnorm(mean=0,sigma=1,upper=y[1]/sig0)[[1]] } for (i in 2:ntot) { cxx <- rho^(jd[i]-jd[i-1])/del[i]/del[i-1] if (!icen[i]) { ztmp <- (y[i]-cxx*yimp[i-1])/(1-cxx^2)^0.5/sig0 if(ztmp>2.5) {ztmp <- 2.5 } else if (ztmp < (-2.5)) {ztmp <- (-2.5) } cumlik <- cumlik + ztmp^2 + 2*log(sig0) + log(1-cxx^2) yimp[i] <- ztmp*(1-cxx^2)^0.5*sig0 + cxx*yimp[i-1] } else { ztmp <- (y[i]-cxx*yimp[i-1])/(1-cxx^2)^0.5/sig0 if(ztmp<(-2.5)) {ztmp <- (-2.5) } cumlik <- cumlik-2*log(pnorm(ztmp)) yimp[i] <- mtmvnorm(mean=0,sigma=1,upper=ztmp)[[1]] yimp[i] <- yimp[i]*(1-cxx^2)^0.5*sig0 + cxx*yimp[i-1] } } cumlik <- cumlik + 2*sum(log(del[!icen]))+ 2*sum(log(omg[!icen])) + 500*(log(0.9) - log(rho))^4 + 0.5*(alph/2)^4 cumlik } ######################################################################################################### compwaveconvxx <- function(cmax,jmod,hlife,spmodcls) { # computes seasonal wave with cmax=time of peak, jmod=model choice (1:6), # and hlife=halflife, in months (1:4) # spmodcls=1 for "one-hump" models, spmodcls=2 for "two-hump" models del <- 1/12 txx <- seq(0,1,1/360) sxx <- txx if (spmodcls==1) { phi <- 12/hlife if(jmod==1) {wtx <- c(1,0,0,0,0,0,0,0,0,0,0,0); pkt <- 1/12} if(jmod==2) {wtx <- c(1,1,0,0,0,0,0,0,0,0,0,0); pkt <- 2/12} if(jmod==3) {wtx <- c(1,1,1,0,0,0,0,0,0,0,0,0); pkt <- 3/12} if(jmod==4) {wtx <- c(1,1,1,1,0,0,0,0,0,0,0,0); pkt <- 4/12} if(jmod==5) {wtx <- c(1,1,1,1,1,1,0,0,0,0,0,0); pkt <- 6/12} if(jmod==6) {wtx <- c(1,1,1,1,1,1,1,1,1,0,0,0); pkt <- 9/12} } else if (spmodcls==2) { if (hlife<=2) { phi <- 12/hlife if(jmod==1) {wtx <- c(1,1,0,0,1,1,1,1,0,0,0,0); pkt <- 8/12} if(jmod==4) {wtx <- c(1,1,1,0,0,1,1,1,0,0,0,0); pkt <- 8/12} if(jmod==3) {wtx <- c(1,1,1,1,0,0,1,1,0,0,0,0); pkt <- 4/12} if(jmod==2) {wtx <- c(1,1,0,0,0,1,1,1,1,0,0,0); pkt <- 9/12} if(jmod==5) {wtx <- c(1,1,1,0,0,0,1,1,1,0,0,0); pkt <- 9/12} if(jmod==6) {wtx <- c(1,1,1,1,0,0,1,1,1,1,0,0); pkt <- 10/12} } else { phi <- 12/(hlife-2) if(jmod==1) {wtx <- c(2,2,0,0,1,1,1,1,0,0,0,0); pkt <- 2/12} if(jmod==4) {wtx <- c(1.5,1.5,1.5,0,0,1,1,1,0,0,0,0); pkt <- 3/12} if(jmod==3) {wtx <- c(1,1,1,1,0,0,2,2,0,0,0,0); pkt <- 8/12} if(jmod==2) {wtx <- c(2,2,0,0,0,1,1,1,1,0,0,0); pkt <- 2/12} if(jmod==5) {wtx <- c(1.5,1.5,1.5,0,0,0,1,1,1,0,0,0); pkt <- 3/12} if(jmod==6) {wtx <- c(1.5,1.5,1.5,1.5,0,0,1,1,1,1,0,0); pkt <- 4/12} } } rho <- exp(-phi) z0xx <- rho^(txx) r12 <- rho^(-del*c(1:12)) con <- wtx[1]*(r12[1]-1) for(k in 2:12) {con <- con+wtx[k]*(r12[k]-r12[k-1])} con <- rho/(1-rho)*con z0xx <- z0xx*con zmat <- matrix(nrow=length(txx),ncol=12) pckm <- matrix(nrow=length(txx),ncol=12) ntot <- length(txx) pckm[,1] <- c(txx<=del) for (k in 2:12) {pckm[,k] <- c(txx>(k-1)*del & txx<=k*del)} zmat[,1] <- rho^txx*(rho^(-replace(txx,txx>1/12,1/12))-1) for (k in 2:12) { ztmp <- rep(0,ntot) for (j in 1:(k-1)) {ztmp[pckm[,j]] <- 0} ztmp[pckm[,k]] <- 1-rho^(txx[pckm[,k]]-(k-1)*del) if(k<12){ for (j in (k+1):12) {ztmp[pckm[,j]] <- rho^(txx[pckm[,j]]-del*k)-rho^(txx[pckm[,j]]-(k-1)*del)} } zmat[,k] <- ztmp } sst <- z0xx for (k in 1:12) { sst <- sst+wtx[k]*zmat[,k]} sst <- sst/phi medxx <- (max(sst)+min(sst))/2 rngxx <- 2*(max(sst)-medxx) sst <- (sst-medxx)/rngxx ###################### if(cmax<=pkt) { txx2 <- (txx+1-pkt+cmax)} if(cmax>pkt) {txx2 <- txx-pkt+cmax} txx2[txx2>1] <- txx2[txx2>1]-1 otmp <- order(txx2); sst <- sst[otmp] sst } ######################################################################################################### swaveqexCSIM <- function(ngen,parxxx,datobs,datall,spmodcls) { # generates conditional simulations from seawaveqex model ## define seasonal wave parameters etc hlife <- parxxx[10] ; cmax <- parxxx[8] ; modn <- parxxx[9] swave <- compwaveconvxx(cmax,modn,hlife,spmodcls) sdum <- (1:365)/366 ipk <- floor(360*sdum) ipk[ipk==0] <- 1 swave <- swave[ipk] waveobs <- swave[datobs[,"sday"]] waveall <- swave[datall[,"sday"]] datobs[,"swave"] <- waveobs datall[,"swave"] <- waveall ## define other model parameters etc. alph <- parxxx[12] ; rho <- parxxx[13] ; sig <- parxxx[14] omgobs <- (1+alph*(0.0+waveobs)) omgall <- (1+alph*(0.0+waveall)) int0 <- parxxx[1] ; wc0 <- parxxx[2]; mtc0 <- parxxx[3]; stc0 <- parxxx[4]; tndc0 <- parxxx[5] muobs <- int0*datobs[,"int"]+wc0*datobs[,"swave"]+mtc0*datobs[,"qanom30"]+stc0*datobs[,"qanom1"]+tndc0*datobs[,"tnd"] muall <- int0*datall[,"int"]+wc0*datall[,"swave"]+mtc0*datall[,"qanom30"]+stc0*datall[,"qanom1"]+tndc0*datall[,"tnd"] yhatobs <- mtc0*datobs[,"qanom30"]+stc0*datobs[,"qanom1"]+tndc0*datobs[,"tnd"] yhatall <- mtc0*datall[,"qanom30"]+stc0*datall[,"qanom1"]+tndc0*datall[,"tnd"] yadj2obs <- wc0*datobs[,"swave"]+tndc0*datobs[,"tnd"]+stc0*datobs[,"qanom1"] yadj3obs <- wc0*datobs[,"swave"]+tndc0*datobs[,"tnd"]+mtc0*datobs[,"qanom30"] yadj2all <- wc0*datall[,"swave"]+tndc0*datall[,"tnd"]+stc0*datall[,"qanom1"] yadj3all <- wc0*datall[,"swave"]+tndc0*datall[,"tnd"]+mtc0*datall[,"qanom30"] yadj4obs <- wc0*datobs[,"swave"]+mtc0*datobs[,"qanom30"]+stc0*datobs[,"qanom1"] yadj4all <- wc0*datall[,"swave"]+mtc0*datall[,"qanom30"]+stc0*datall[,"qanom1"] yobs <- datobs[,"yobs"] ; ytilobs <- (yobs-muobs)/omgobs yall <- datall[,"yobs"] ; ytilall <- (yall-muall)/omgall jdayobs <- datobs[,"jday"] jdayall <- datall[,"jday"] icen <- datobs[,"icen"]=="c" svcondsim <- matrix(nrow=length(jdayall),ncol=ngen) for (igen in 1:ngen) { # impute censored values if necessary yimp <- ytilobs if (sum(icen) > 0) { yimp <- impcenvals(rho,sig,yimp,icen,jdayobs) } # now do conditional simulation ysim <- condsim(rho,sig,yimp,jdayobs,jdayall) ysim <- ysim*omgall + muall yimp <- yimp*omgobs + muobs svcondsim[,igen] <- ysim } ## compute various statistics etc for output ... svcondmean <- apply(10^svcondsim,1,mean) svcondmean <- log10(svcondmean) svcondsim1 <- svcondsim[,1] nall <- length(jdayall) nyrs <- round(nall/365,0) tmpannmax <- matrix(nrow=nyrs,ncol=ngen) flagnas <- rep(F,nyrs) for (i in 1:ngen) { tmpyy <- 10^svcondsim[,i] tmpyy <- matrix(tmpyy,nrow=365) for (j in 1:nyrs) { if (sum(!is.na(tmpyy[,j]))<90) { tmpannmax[j,i] <- NA flagnas[j] <- T } else if(sum(!is.na(tmpyy[,j]))<365) { tmpannmax[j,i] <- max(tmpyy[,j],na.rm=T) flagnas[j] <- T } else { tmpannmax[j,i] <- max(tmpyy[,j]) } } } maxout <- apply(tmpannmax,1,mean) #### **** save 10th and 90th percentiles p10out <- apply(tmpannmax,1,quantile,probs=0.1,na.rm=T) p90out <- apply(tmpannmax,1,quantile,probs=0.9,na.rm=T) p10out[is.na(maxout)] <- NA p90out[is.na(maxout)] <- NA modoutobs <- data.frame(yhatobs,muobs,omgobs,svcondsim1[jdayobs],yadj2obs,yadj3obs,yadj4obs) names(modoutobs) <- c("yhat","mu","omg","condsim","yadj2","yadj3","yadj4") modoutall <- data.frame(yhatall,muall,omgall,svcondsim1,svcondmean,yadj2all,yadj3all,yadj4all) names(modoutall) <- c("yhat","mu","omg","condsim","condmean","yadj2","yadj3","yadj4") list(maxout,modoutobs,modoutall,svcondsim,p10out,p90out,flagnas) } ######################################################################################################### impcenvals <- function(rho,sig,yimp,icen,jdayobs) { # # generates conditional values for censored data only # nobs <- length(yimp) ncum <- 1:nobs indu <- !icen ; iu <- ncum[indu] nu <- sum(indu) nblkc <- nu+1 ncb <- rep(0,nblkc) if(icen[1]) {ncb[1] <- iu[1]-1} for (i in 2:nu) {ncb[i] <- iu[i]-iu[i-1]-1} if (icen[nobs]) {ncb[nu+1] <- nobs-iu[nu]} t1 <- jdayobs #### ready to impute new values for censored data for (i in 1:(nu+1)) { if(i==1 & ncb[1]>0) { yblk <- yimp[1:ncb[1]] yends <- yimp[ncb[1]+1] mtmp <- matrix(nrow=ncb[1]+1,ncol=ncb[1]+1) for (m in 1:(ncb[1]+1)) { for (k in m:(ncb[1]+1)) { if(k==m) {mtmp[k,k] <- 1} if(k>m) { mtmp[m,k] <- rho^(t1[k]-t1[m]) mtmp[k,m] <- mtmp[m,k] } } } mtmpx <- as.matrix(mtmp[1:ncb[1],1:ncb[1]]); gtmpx <- as.matrix(mtmp[1:ncb[1],ncb[1]+1],ncol=1) mtmpy <- mtmp[ncb[1]+1,ncb[1]+1] cvar <- mtmpx-gtmpx%*%t(gtmpx)/mtmpy # cvarinv <- solve(cvar) cmu <- gtmpx%*%as.matrix(yends,ncol=1)/mtmpy yimp[1:ncb[1]] <- rtmvnorm(n=1,mean=c(cmu),,sigma=sig^2*cvar,upper=yblk,algorithm='gibbs',burn.in.samples=100) } if (i>1 & i<(nu+1) & ncb[i]>0) { yblk <- yimp[(iu[i-1]+1):(iu[i]-1)] yends <- yimp[c(iu[i-1],iu[i])] mtmp <- matrix(nrow=ncb[i]+2,ncol=ncb[i]+2) for (m in 1:(ncb[i])) { for (k in m:(ncb[i])) { if(k==m) {mtmp[k,k] <- 1} if(k>m) { mtmp[m,k] <- rho^(t1[iu[i-1]+k]-t1[iu[i-1]+m]) mtmp[k,m] <- mtmp[m,k] } } } for (m in 1:ncb[i]) { mtmp[m,ncb[i]+1] <- rho^(t1[iu[i-1]+m]-t1[iu[i-1]]) mtmp[ncb[i]+1,m] <- mtmp[m,ncb[i]+1] mtmp[m,ncb[i]+2] <- rho^(t1[iu[i]]-t1[iu[i-1]+m]) mtmp[ncb[i]+2,m] <- mtmp[m,ncb[i]+2] } mtmp[ncb[i]+1,ncb[i]+1] <- 1 mtmp[ncb[i]+2,ncb[i]+2] <- 1 mtmp[ncb[i]+1,ncb[i]+2] <- rho^(t1[iu[i]]-t1[iu[i-1]]) mtmp[ncb[i]+2,ncb[i]+1] <- mtmp[ncb[i]+1,ncb[i]+2] mtmpx <- as.matrix(mtmp[1:ncb[i],1:ncb[i]]); gtmpx <- matrix(mtmp[1:ncb[i],c(ncb[i]+1,ncb[i]+2)],ncol=2) mtmpy <- mtmp[c(ncb[i]+1,ncb[i]+2),c(ncb[i]+1,ncb[i]+2)] cvar <- mtmpx-gtmpx%*%solve(mtmpy)%*%t(gtmpx) cmu <- gtmpx%*%solve(mtmpy)%*%as.matrix(yends,ncol=1) # cvarinv <- solve(cvar) yimp[(iu[i-1]+1):(iu[i-1]+ncb[i])] <- rtmvnorm(n=1,mean=c(cmu),sigma=sig^2*cvar,upper=yblk,algorithm='gibbs',burn.in.samples=100) } if(i==(nu+1) & ncb[nu+1]>0) { yblk <- yimp[(iu[nu]+1):(iu[nu]+ncb[nu+1])] yends <- yimp[iu[nu]] mtmp <- matrix(nrow=ncb[nu+1]+1,ncol=ncb[nu+1]+1) for (m in 1:(ncb[nu+1]+1)) { for (k in m:(ncb[nu+1]+1)) { if(k==m) {mtmp[k,k] <- 1} if(k>m & k<(ncb[nu+1]+1)) { mtmp[m,k] <- rho^(t1[iu[nu]+k]-t1[iu[nu]+m]) mtmp[k,m] <- mtmp[m,k] } if(k>m & k==(ncb[nu+1]+1)) { mtmp[m,k] <- rho^(t1[iu[nu]+m]-t1[iu[nu]]) mtmp[k,m] <- mtmp[m,k] } } } mtmpx <- as.matrix(mtmp[1:ncb[nu+1],1:ncb[nu+1]]); gtmpx <- as.matrix(mtmp[1:ncb[nu+1],ncb[nu+1]+1],ncol=1) mtmpy <- mtmp[ncb[nu+1]+1,ncb[nu+1]+1] cvar <- mtmpx-gtmpx%*%t(gtmpx)/mtmpy cmu <- gtmpx%*%as.matrix(yends,nrow=1)/mtmpy # cvarinv <- solve(cvar) yimp[(iu[nu]+1):(iu[nu]+ncb[nu+1])] <- rtmvnorm(n=1,mean=c(cmu),sigma=sig^2*cvar,upper=yblk,algorithm='gibbs',burn.in.samples=100) } } yimp[yimp<(-3)] <- (-3) yimp } ######################################################################################################### condsim <- function(rho,sig,yimp,jdayobs,jdayall) { # # generates conditional simulations # nobs <- length(yimp) nall <- length(jdayall) ysim <- rep(NA,nall) tobs <- jdayobs tall <- jdayall for (i in 1:nobs) { ysim[tobs[i]] <- yimp[i] if(i==1 & tobs[1]>1) { n <- tobs[1] nm1 <- n-1 yends <- yimp[1] mtmp <- matrix(nrow=n,ncol=n) for (m in 1:n) { for (k in m:n) { if(k==m) {mtmp[k,k] <- 1} if(k>m) { mtmp[m,k] <- rho^(k-m) mtmp[k,m] <- mtmp[m,k] } } } mtmpx <- as.matrix(mtmp[1:nm1,1:nm1]); gtmpx <- as.matrix(mtmp[1:nm1,n]) mtmpy <- mtmp[n,n] cvar <- mtmpx-gtmpx%*%t(gtmpx)/mtmpy # cvarinv <- solve(cvar) cmu <- gtmpx%*%as.matrix(yends,ncol=1)/mtmpy ysim[1:nm1] <- rmvnorm(n=1,mean=cmu,sigma=sig^2*cvar,method="chol") } if (i>1 & i<=nobs & (tobs[i]-tobs[max(1,i-1)])>0) { n <- tobs[i]-tobs[i-1]+1 nm2 <- n-2 yends <- yimp[c(i-1,i)] mtmp <- matrix(nrow=n,ncol=n) for (m in 1:nm2) { for (k in m:nm2) { if(k==m) {mtmp[k,k] <- 1} if(k>m) { mtmp[m,k] <- rho^(k-m) mtmp[k,m] <- mtmp[m,k] } } } for (m in 1:nm2) { mtmp[m,nm2+1] <- rho^(m) mtmp[nm2+1,m] <- mtmp[m,nm2+1] mtmp[m,nm2+2] <- rho^(nm2-m+1) mtmp[nm2+2,m] <- mtmp[m,nm2+2] } mtmp[nm2+1,nm2+1] <- 1 mtmp[nm2+2,nm2+2] <- 1 mtmp[nm2+1,nm2+2] <- rho^(tobs[i]-tobs[i-1]) mtmp[nm2+2,nm2+1] <- mtmp[nm2+1,nm2+2] mtmpx <- as.matrix(mtmp[1:nm2,1:nm2]); gtmpx <- matrix(mtmp[1:nm2,c(nm2+1,nm2+2)],ncol=2) mtmpy <- mtmp[c(nm2+1,nm2+2),c(nm2+1,nm2+2)] cvar <- mtmpx-gtmpx%*%solve(mtmpy)%*%t(gtmpx) cmu <- gtmpx%*%solve(mtmpy)%*%as.matrix(yends,ncol=1) ysim[(tobs[i-1]+1):(tobs[i]-1)] <- rmvnorm(n=1,mean=c(cmu),sigma=sig^2*cvar,method="chol") } if(i==nobs & tobs[i]m) { mtmp[m,k] <- rho^(k-m) mtmp[k,m] <- mtmp[m,k] } } } for (m in 1:nm1) { mtmp[m,n] <- rho^(m) mtmp[n,m] <- mtmp[m,n] } mtmp[n,n] <- 1 mtmpx <- as.matrix(mtmp[1:nm1,1:nm1]); gtmpx <- as.matrix(mtmp[1:nm1,n]) mtmpy <- mtmp[n,n] cvar <- mtmpx-gtmpx%*%t(gtmpx)/mtmpy # cvarinv <- solve(cvar) cmu <- gtmpx%*%as.matrix(yends,ncol=1)/mtmpy ysim[(tobs[i]+1):nall] <- rmvnorm(n=1,mean=cmu,sigma=sig^2*cvar,method="chol") } } ysim[ysim<(-3)] <- (-3) ysim } ######################################################################## that's it!