m <- as.matrix(CIframe.by.age.angio.plus.recent3)

cor <- rep (0, 15)
rho <- rep (0, 15)
tau <- rep (0, 15)
#p <- rep (0, 15)
#r2 <- rep (0, 15)

corcilow <- rep (0, 15)
corcihi <- rep (0, 15)

Rcor <- rep (0, 16)
Rrho <- rep (0, 16)
Rtau <- rep (0, 16)
#Rp <- rep (0, 16)
#Rr2 <- rep (0, 16)

Rcorcilow <- rep (0, 16)
Rcorcihi <- rep (0, 16)

for(i in 2:16){
 #m.lm <- lm (m[i, ] ~ m[i-1, ])
 cor[i-1] <- format(cor.test (as.numeric(m[i, ]), as.numeric(m[(i-1), ]),
		method = 'pearson', alternative = 'greater')$estimate, digits = 2)
 corcilow[i-1] <- as.numeric(format(cor.test (as.numeric(m[i, ]),
		as.numeric(m[(i-1), ]), method = 'pearson',
		alternative = 'greater')$conf.int[1], digits = 2))
 corcihi[i-1] <- as.numeric(format(cor.test (as.numeric(m[i, ]),
		as.numeric(m[(i-1), ]), method = 'pearson',
		alternative = 'greater')$conf.int[2], digits = 2))

 rho[i-1] <- format (cor.test (m[i, ], m[(i-1), ], method = 'spearman',
		alternative = 'greater')$estimate, digits = 2)
 tau[i-1] <- format (cor.test (m[i, ], m[(i-1), ], method = 'kendall',
		alternative = 'greater')$estimate, digits = 2)
 #Rm.lm <- lm (m[1, ] ~ m[i-1, ])
 Rcor[i-1] <- format (cor.test (as.numeric(m[1, ]), as.numeric(m[(i-1), ]),
		method = 'pearson', alternative = 'greater')$estimate, digits = 2)
 Rcorcilow[i-1] <- cor.test (as.numeric(m[1, ]), as.numeric(m[(i-1), ]),
		method = 'pearson', alternative = 'greater')$conf.int[1]
 Rcorcihi[i-1] <- cor.test (as.numeric(m[1, ]), as.numeric(m[(i-1), ]),
		method = 'pearson', alternative = 'greater')$conf.int[2]

 Rrho[i-1] <- format (cor.test (m[1, ], m[(i-1), ], method = 'spearman',
		alternative = 'greater')$estimate, digits = 2)
 Rtau[i-1] <- format (cor.test (m[1, ], m[(i-1), ], method = 'kendall',
		alternative = 'greater')$estimate, digits = 2)
}

Rcor[i] <- format (cor.test (as.numeric(m[1, ]), as.numeric(m[i, ]),
		method = 'pearson', alternative = 'greater')$estimate, digits = 2)
#Rcorcilow[i] <- cor.test (as.numeric(m[1, ]), as.numeric(m[i, ]),
		method = 'pearson', alternative = 'greater')$conf.int[1]
#Rcorcihi[i] <- cor.test (as.numeric(m[1, ]), as.numeric(m[i, ]),
		method = 'pearson', alternative = 'greater')$conf.int[2]
Rcorcilow[i] <- -1
Rcorcihi[i] <- 1
Rrho[i] <- format (cor.test (m[1, ], m[i, ], method = 'spearman',
		alternative = 'greater')$estimate, digits = 2)
Rtau[i] <- format (cor.test (m[1, ], m[i, ], method = 'kendall',
		alternative = 'greater')$estimate, digits = 2)

# This has produced 15 autocorrelations (for the 16 time intervals) and
# 16 correlations to the recent (incuding the trivial first one, which is equal to 1).

boundaries <- c(0, ChronoStratDates$GSAStartDate[1:16]) # 17 boundaries for Barremian--present
centers <- c (mean(boundaries[1:2]), mean(boundaries[2:3]),
		mean(boundaries[3:4]), mean(boundaries[4:5]),
		mean(boundaries[5:6]), mean(boundaries[6:7]),
		mean(boundaries[7:8]), mean(boundaries[8:9]),
		mean(boundaries[9:10]), mean(boundaries[10:11]),
		mean(boundaries[11:12]), mean(boundaries[12:13]),
		mean(boundaries[13:14]), mean(boundaries[14:15]),
		mean(boundaries[15:16]), mean(boundaries[16:17])) # 16 centers for these time intervals

#boundaries
#centers

#library(vegan)
#horndist <- vegdist(m, method = 'horn')

pdf(file = 'CorCoefThroughTime.pdf', width = 8, height = 6, family = 'Times', pointsize = 10)
par(mfrow = c(2,1))

plot.default (-(boundaries[2:16]), cor[1:15], type = 'l', lty = 1, xaxt = 'n', ann = FALSE, xlim = c(-124, 0), ylim = c(-0.2, 1))

title(ylab = 'Correlation', main = 'Autocorrelation')

lines (-(boundaries[2:16]), corcilow[1:15], lty = 3, lwd = 0.3)
lines (-(boundaries[2:16]), corcihi[1:15], lty = 3, lwd = 0.3)

lines (-(boundaries[2:16]), rho[1:15], col = 2, lty = 2)

lines (-(boundaries[2:16]), tau[1:15], col = 3, lty = 4)

#lines (as.matrix(horndist)[16,], type = 'l', col = 'red')

axis (1, at = -(centers[1:16]), labels = row.names(m)[1:16], tick = FALSE)
axis (1, at = -(boundaries[2:16]), labels = FALSE, tick = TRUE)
#axis (1, at = 0.5:15.5, labels = row.names(m)[1:16])

abline (v = -(boundaries[8]), col = 'red', lwd = 3)

corbar <- rowSums(cbind (as.numeric(cor), as.numeric(rho), as.numeric(tau)))/3
abline (lm (corbar ~ boundaries[2:16])$coef[1], -lm(corbar ~ boundaries[2:16])$coef[2], col = 7)

summary(lm (corbar ~ boundaries[2:16]))
mtext (paste ('trendp', summary(lm (corbar ~ boundaries[2:16]))$coef[2,4], sep = ' = '))

plot.default (-(centers), Rcor[1:16], type = 'n', xaxt = 'n', ann = FALSE, xlim = c(-124, 0), ylim = c(-0.2, 1))

title(ylab = 'Correlation', main = 'Correlation to Present')

lines (-(centers), Rcorcilow[1:16], lty = 3, lwd = 0.3)
lines (-(centers), Rcorcihi[1:16], lty = 3, lwd = 0.3)

lines (-(centers), Rcor[1:16], col = 1, lty = 1)
lines (-(centers), Rrho[1:16], col = 2, lty = 2)
lines (-(centers), Rtau[1:16], col = 3, lty = 4)

#lines (as.matrix(horndist)[,1], type = 'l', col = 'red')

axis (1, at = -(centers[1:16]), labels = row.names(m)[1:16], tick = FALSE)
axis (1, at = -(boundaries[2:16]), labels = FALSE, tick = TRUE)

abline (v = -(boundaries[8]), col = 'red', lwd = 3)

Rcorbar <- rowSums(cbind (as.numeric(Rcor), as.numeric(Rrho), as.numeric(Rtau)))/3
abline (lm (Rcorbar ~ centers)$coef[1], -lm(Rcorbar ~ centers)$coef[2], col = 7)
abline (lm (Rcorbar[2:14] ~ centers[2:14])$coef[1], -lm(Rcorbar[2:14] ~ centers[2:14])$coef[2], col = 6)

summary(lm (Rcorbar ~ centers))
mtext (paste ('trendp', summary(lm (Rcorbar ~ centers))$coef[2,4], sep = ' = '))
mtext (paste ('trendpminusoutliers', summary(lm (Rcorbar[2:14] ~ centers[2:14]))$coef[2,4], sep = ' = '))

dev.off()


#macintosh()
#par (mfcol = c(2,2))
#plot(lm (corbar ~ boundaries[2:16]))
#dev.off()

#macintosh()
#par (mfcol = c(2,2))
#plot(lm (corbar[1:14] ~ boundaries[2:15]))
#dev.off()

#macintosh()
#par (mfcol = c(2,2))
#plot(lm (Rcorbar ~ centers))
#dev.off()

#macintosh()
#par (mfcol = c(2,2))
#plot(lm (Rcorbar[2:14] ~ centers[2:14]))
#dev.off()

pdf(file = 'CorCoefDists.pdf', width = 5, height = 5, family = 'Times')
#macintosh()
par (mfcol = c(3,1))
hist (as.numeric(cor), xlim = c(-0.2, 1), main = 'Distribution of Correlation Coefficients Between Adjacent Ages', xlab = 'Pearson\'s Product Momement Correlation Coefficient', breaks = seq (-0.4, 1, by = 0.1))
lines (density (as.numeric(cor)))
abline (v = as.numeric(cor[9]), col = 2, lwd = 2)
hist (as.numeric(rho), border = 2, xlim = c(-0.2, 1), main = 'Distribution of Correlation Coefficients Between Adjacent Ages', xlab = 'Spearman\'s Rank-Order Correlation Coefficient Rho', breaks = seq (-0.4, 1, by = 0.1))
lines (density (as.numeric(rho)), col = 2)
abline (v = as.numeric(rho[9]), col = 2, lwd = 2)
hist (as.numeric(tau), border = 3, xlim = c(-0.2, 1), main = 'Distribution of Correlation Coefficients Between Adjacent Ages', xlab = 'Kendall\'s Rank-Order Correlation Coefficient Tau', breaks = seq (-0.4, 1, by = 0.1))
lines (density (as.numeric(tau)), col = 3)
abline (v = as.numeric(tau[9]), col = 2, lwd = 2)
dev.off()

pdf(file = 'MaasPaleoBiplot.pdf', width = 4, height = 4, pointsize = 11, family = 'Times')
p <- CI$frame.by.age[1:56,6]
m <- CI$frame.by.age[1:56,7]

p[is.na(p)] <- 0
m[is.na(m)] <- 0

vegdist(rbind(m,p), method = 'horn')

p[p == 0 & m == 0] <- NA
m[is.na(p) & m == 0] <- NA

summary (lm(p ~ m))

plot.default (m, p, type = 'n', asp = 1, xlab = 'Maastrictian Counts', ylab = 'Paleocene Counts', main = 'Comparison of counts in each CIC')
text (m, p, labels = CICs[1:56])
abline (coef = c(0, 1))
abline (lm(p ~ m)$coef, lty = 2)
mtext (paste ('p', summary(lm (p ~ m))$coef[2,4], sep = ' = '))
rm (m, p)
dev.off()

#macintosh()
par (mfcol = c(2,2))
plot(lm (p ~ m))
plot (lm((p[c(-4, -5, -7)]) ~ (m[c(-4, -5, -7)])), cex.id = 1.5, labels.id = CICs[c(2,3,4,5,7, 9, 12, 13)])
plot (lm((p) ~ (m)), cex.id = 1.5, labels.id = CICs[c(2,3,4,5,7, 9, 12, 13)])
plot (lm(log(p) ~ log(m)), cex.id = 1.5, labels.id = CICs[c(2,3,4,5,7, 9, 12, 13)])
#dev.off()

summary (lm((p[c(-4, -5, -7)]) ~ (m[c(-4, -5, -7)])))

#Call:
#lm(formula = (p[c(-4, -5, -7)]) ~ (m[c(-4, -5, -7)]))

#Residuals:
#      2       3       6       9      10 
# 1.5353  0.1235 -0.7294 -1.2000  0.2706 

#Coefficients:
#                 Estimate Std. Error t value Pr(>|t|)  
#(Intercept)        0.9941     0.8272   1.202   0.3157  
#m[c(-4, -5, -7)]   0.7353     0.2081   3.533   0.0386 *
#---
#Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

#Residual standard error: 1.213 on 3 degrees of freedom
#Multiple R-Squared: 0.8062,	Adjusted R-squared: 0.7417 
#F-statistic: 12.48 on 1 and 3 DF,  p-value: 0.03855 

# can also do the above with logged m and p; results are similar and heteroscedasticity is avoided

cor.test(as.numeric(CIframe.by.age.angio.plus.recent3[7,]), as.numeric(CIframe.by.age.angio.plus.recent3[8,]), alternative = 'greater', method = 'pearson')

cor.test(as.numeric(CIframe.by.age.angio.plus.recent3[7,]), as.numeric(CIframe.by.age.angio.plus.recent3[8,]), alternative = 'greater', method = 'spearman)

cor.test(as.numeric(CIframe.by.age.angio.plus.recent3[7,]), as.numeric(CIframe.by.age.angio.plus.recent3[8,]), alternative = 'greater', method = 'kendall')

rm (m, i, cor, rho, tau, p, r2, m.lm, Rcor, Rrho, Rtau, Rp, Rr2, Rm.lm, corcilow, corcihi, Rcorcilow, Rcorcihi, boundaries, centers, corbar, Rcorbar, p)