Data, replicated analyses/plots and other trend tests in r

Reckhowetal_trend.r

## Data from Reckhow, K.H., K Kepford and W.W. Hicks. 1993
### Methods for the Analysis of Lake Water Quality
### EPA 841-R-93-003

## Code was compiled by Paul Julian
## contact info: paul.julian@floridadep.gov



# Appendix C --------------------------------------------------------------

fall.lake<-data.frame(
  OBS=1:60,
  DATE=as.Date(c("1983-01-26","1983-02-25","1983-03-27","1983-04-26","1983-05-19","1983-06-28","1983-07-28","1983-08-31","1983-09-27","1983-10-26","1983-11-29","1983-12-14",
                 "1984-01-09","1984-02-09","1984-03-14","1984-04-17","1984-05-15","1984-06-14","1984-07-19","1984-08-16","1984-09-17","1984-10-23","1984-11-15","1984-12-18",
                 "1985-01-30","1985-02-21","1985-03-20","1985-04-17","1985-05-01","1985-06-12","1985-07-10","1985-08-13","1985-09-24","1985-10-23","1985-11-14","1985-12-27",
                 "1986-01-16","1986-02-20","1986-03-12","1986-04-28","1986-05-27","1986-06-24","1986-07-16","1986-08-06","1986-09-30","1986-10-28","1986-11-18","1986-12-15",
                 "1987-01-15","1987-02-03","1987-03-25","1987-04-07","1987-05-13","1987-06-04","1987-07-01","1987-08-12","1987-09-15","1987-10-14","1987-11-13","1987-12-13")),
  TP=c(NA,NA,NA,0.07,0.04,0.05,0.04,0.04,0.03,0.03,0.03,NA,0.07,0.09,0.10,0.07,0.04,0.03,0.02,0.05,0.04,0.04,0.04,0.03,0.03,0.09,0.05,0.03,0.03,0.02,0.04,0.03,0.02,0.02,0.03,0.05,0.04,0.06,0.06,0.04,0.04,0.03,0.02,0.03,0.01,0.02,0.03,0.02,NA,0.08,0.08,0.04,0.04,0.03,0.04,0.02,0.03,0.02,NA,NA),
  TN=c(NA,NA,NA,0.58,0.36,0.36,0.31,0.51,0.42,0.71,0.87,NA,0.85,0.76,0.72,0.65,0.45,0.41,0.41,0.41,0.51,0.62,0.73,0.63,0.85,0.93,0.71,0.50,0.33,0.41,0.41,0.51,0.41,0.73,0.94,0.74,0.76,0.87,0.80,0.40,0.31,0.41,0.41,0.41,0.41,0.61,0.94,0.88,NA,0.85,0.55,0.44,0.41,0.41,0.21,0.41,0.52,0.81,NA,NA)
)
fall.lake$month<-format(fall.lake$DATE,"%m")
fall.lake$year<-format(fall.lake$DATE,"%Y")



# Results -----------------------------------------------------------------

## Page 23
library(e1071)
# Moments
length(which(fall.lake$TP!="NA")); #N
mean(fall.lake$TP,na.rm=T); #mean
sd(fall.lake$TP,na.rm=T); #standard deviation
e1071::skewness(fall.lake$TP,na.rm=T,type=2); #skewness
sum(fall.lake$TP^2,na.rm=TRUE); #corrected sum of squares (USS)
(sd(fall.lake$TP,na.rm=T)/mean(fall.lake$TP,na.rm=T))*100; #coefficient of variance (CV)
sum(fall.lake$TP,na.rm=T)
var(fall.lake$TP,na.rm=T)
e1071::kurtosis(fall.lake$TP,na.rm=T,type=2)
sum((fall.lake$TP-mean(fall.lake$TP,na.rm=T))^2,na.rm=TRUE);#corrected sum of squares (CSS)

# omitted T, S and W statistics. Not clear what  these are.
#quantiles
round(quantile(fall.lake$TP,na.rm=T,probs=c(0,0.01,0.05,0.10,0.25,0.50,0.75,0.90,0.95,0.99,1)),2)

diff(range(fall.lake$TP,na.rm=T)); #range
diff(quantile(fall.lake$TP,na.rm=T,probs=c(0.25,0.75))); #inter-quantile range

## Page 24
boxplot(fall.lake$TP)

qqnorm(fall.lake$TP)
qqline(fall.lake$TP,col="red")

## Page 26
hist(fall.lake$TP,xlim=c(0,0.1))

## Page 27
plot(TP~DATE,fall.lake,type="b")
abline(lm(TP~DATE,fall.lake),lty=3); #despite this being incorrect

## Page 28
boxplot(TP~month,fall.lake,notch=T,ylim=c(0,0.12))

## Page 30
boxplot(TP~year,fall.lake,notch=T,ylim=c(0,0.10))

## Page 32
acf(na.omit(fall.lake$TP))
acf(na.omit(fall.lake$TP),plot=F)


########
########
# Different Trend Tests ---------------------------------------------------
########
########
## Page 36 of report (Table 3.4)
#tau statistic -0.29787
#p-value without serial correlation 0.013974
#p-value with serial correlation 0.23909
#slope statistic -0.0033333

# Base R for standard kendall
cor.test(fall.lake$OBS,fall.lake$TP,method="kendall")

# as a ts.object
fall.lake.ts<-ts(fall.lake$TP,frequency=12,start=c(1983,1))

plot(fall.lake.ts)



# Regional/Seasonal Kendall consistent with USGS methods
library(rkt)
library(lubridate)
with(na.omit(fall.lake),rkt(decimal_date(DATE),TP,as.numeric(month)))

# Thiel Sen estimators
library(mblm)
mblm(TP~OBS,na.omit(fall.lake),repeated = F)

library(zyp)
zyp.sen(TP~OBS,na.omit(fall.lake))

library(Kendall)
summary(with(na.omit(fall.lake),Kendall(OBS,TP)))

SeasonalMannKendall(fall.lake.ts)