3
$\begingroup$

I'm struggling with figuring out the best way to break up my river discharge data in a way that I can use it as a factor for further analysis. I'm currently manipulating the data in R but plan to pull it into PRIMER to create a multi-dimensional scaling plot and other analysis.

I have the following data that breaks down by "Zone" (River B and River D), "Year", "Month", and "Discharge" (cubic ft/sec).

structure(list(ï..Zone = c("D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", 
"D", "D", "D", "D", "D", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", 
"B", "B", "B", "B"), Year = c(2005L, 2006L, 2007L, 2008L, 2009L, 
2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 2018L, 
2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 
2014L, 2015L, 2016L, 2017L, 2018L, 2005L, 2006L, 2007L, 2008L, 
2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 
2018L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 
2013L, 2014L, 2015L, 2016L, 2017L, 2018L, 2005L, 2006L, 2007L, 
2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 
2017L, 2018L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 
2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 2018L, 2005L, 2006L, 
2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 
2016L, 2017L, 2018L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 
2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 2018L, 2005L, 
2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 
2015L, 2016L, 2017L, 2018L, 2005L, 2006L, 2007L, 2008L, 2009L, 
2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 2018L, 
2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 
2014L, 2015L, 2016L, 2017L, 2018L, 2005L, 2006L, 2007L, 2008L, 
2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 
2018L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 
2013L, 2014L, 2015L, 2016L, 2017L, 2018L, 2005L, 2006L, 2007L, 
2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 
2017L, 2018L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 
2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 2018L, 2005L, 2006L, 
2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 
2016L, 2017L, 2018L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 
2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 2018L, 2005L, 
2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 
2015L, 2016L, 2017L, 2018L, 2005L, 2006L, 2007L, 2008L, 2009L, 
2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 2018L, 
2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 
2014L, 2015L, 2016L, 2017L, 2018L, 2005L, 2006L, 2007L, 2008L, 
2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 
2018L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 
2013L, 2014L, 2015L, 2016L, 2017L, 2018L, 2005L, 2006L, 2007L, 
2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 
2017L, 2018L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 
2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 2018L), Month = c("Jan", 
"Jan", "Jan", "Jan", "Jan", "Jan", "Jan", "Jan", "Jan", "Jan", 
"Jan", "Jan", "Jan", "Jan", "Feb", "Feb", "Feb", "Feb", "Feb", 
"Feb", "Feb", "Feb", "Feb", "Feb", "Feb", "Feb", "Feb", "Feb", 
"Mar", "Mar", "Mar", "Mar", "Mar", "Mar", "Mar", "Mar", "Mar", 
"Mar", "Mar", "Mar", "Mar", "Mar", "April", "April", "April", 
"April", "April", "April", "April", "April", "April", "April", 
"April", "April", "April", "April", "May", "May", "May", "May", 
"May", "May", "May", "May", "May", "May", "May", "May", "May", 
"May", "June", "June", "June", "June", "June", "June", "June", 
"June", "June", "June", "June", "June", "June", "June", "July", 
"July", "July", "July", "July", "July", "July", "July", "July", 
"July", "July", "July", "July", "July", "August", "August", "August", 
"August", "August", "August", "August", "August", "August", "August", 
"August", "August", "August", "August", "September", "September", 
"September", "September", "September", "September", "September", 
"September", "September", "September", "September", "September", 
"September", "September", "October", "October", "October", "October", 
"October", "October", "October", "October", "October", "October", 
"October", "October", "October", "October", "November", "November", 
"November", "November", "November", "November", "November", "November", 
"November", "November", "November", "November", "November", "November", 
"December", "December", "December", "December", "December", "December", 
"December", "December", "December", "December", "December", "December", 
"December", "December", "Jan", "Jan", "Jan", "Jan", "Jan", "Jan", 
"Jan", "Jan", "Jan", "Jan", "Jan", "Jan", "Jan", "Jan", "Feb", 
"Feb", "Feb", "Feb", "Feb", "Feb", "Feb", "Feb", "Feb", "Feb", 
"Feb", "Feb", "Feb", "Feb", "Mar", "Mar", "Mar", "Mar", "Mar", 
"Mar", "Mar", "Mar", "Mar", "Mar", "Mar", "Mar", "Mar", "Mar", 
"April", "April", "April", "April", "April", "April", "April", 
"April", "April", "April", "April", "April", "April", "April", 
"May", "May", "May", "May", "May", "May", "May", "May", "May", 
"May", "May", "May", "May", "May", "June", "June", "June", "June", 
"June", "June", "June", "June", "June", "June", "June", "June", 
"June", "June", "July", "July", "July", "July", "July", "July", 
"July", "July", "July", "July", "July", "July", "July", "July", 
"August", "August", "August", "August", "August", "August", "August", 
"August", "August", "August", "August", "August", "August", "August", 
"September", "September", "September", "September", "September", 
"September", "September", "September", "September", "September", 
"September", "September", "September", "September", "October", 
"October", "October", "October", "October", "October", "October", 
"October", "October", "October", "October", "October", "October", 
"October", "November", "November", "November", "November", "November", 
"November", "November", "November", "November", "November", "November", 
"November", "November", "November", "December", "December", "December", 
"December", "December", "December", "December", "December", "December", 
"December", "December", "December", "December", "December"), 
    Discharge = c(205.5, 290.1, 58.6, 7.66, 18.5, 465.6, 141.3, 
    3.63, 51.8, 311.8, 210.2, 142.7, 134.7, 133.9, 153.3, 258.6, 
    95.8, 12.7, 26.8, 879.9, 417.3, 3.5, 78.9, 590.6, 347.7, 
    223.3, 75.8, 305.5, 509.3, 99.4, 58.5, 99.7, 23.1, 451.4, 
    167.4, 3.94, 228.7, 1419, 249.9, 148.5, 33.3, 467.2, 1094, 
    26.3, 15.2, 19.2, 164.5, 125.9, 77, 2.59, 173, 1309, 172.3, 
    210, 417.7, 707.4, 489.3, 9.15, 7.88, 4.29, 30.6, 488.5, 
    33.7, 1.88, 52.9, 1191, 53.9, 24.3, 36.8, 278.7, 421.1, 9.94, 
    9.1, 2.47, 59.7, 253.5, 5.25, 620.1, 104.6, 468.1, 32.3, 
    27.5, 592.3, 354.3, 1089, 5.69, 11.6, 5.59, 399.8, 532, 4.31, 
    1508, 1472, 120.6, 299.2, 44.4, 268.1, 849.2, 352.4, 9.98, 
    16.1, 423.2, 122, 382.4, 6.47, 1637, 1545, 86.7, 1099, 1326, 
    535.5, 1066, 40.5, 7.35, 21, 176.3, 120.8, 166.8, 5.67, 1188, 
    724.2, 406, 274.6, 1557, 354.6, 674.4, 16.9, 4.16, 14.5, 
    30.1, 34.1, 57.7, 3.59, 397.4, 164.1, 191.4, 151.1, 108.2, 
    102.8, 153.7, 7.72, 4.29, 8.47, 11.3, 13.3, 14.5, 3.41, 54.3, 
    53.8, 37.5, 54.7, 14.7, 48.9, 141.4, 101.5, 8.95, 8.23, 18.4, 
    92.5, 14.1, 4, 40.3, 67.4, 83.1, 38.5, 27.2, 82.6, 1415, 
    160.3, 286.5, 108.8, 0.13, 31.7, 313, 46.5, 4.46, 40, 261.5, 
    182, 115.5, 44.7, 31.2, 100.5, 280.7, 149.9, 7.89, 34.7, 
    498.2, 123.3, 4.77, 49.4, 283.9, 446.9, 164.2, 39, 50.1, 
    332.8, 222.7, 53.5, 49.5, 29.1, 461.5, 68, 33.8, 189.2, 556.9, 
    260.7, 106.6, 25, 66.4, 684.7, 49.7, 22.7, 28.7, 465.5, 150.8, 
    82.4, 11.7, 198.5, 658.2, 193.5, 260.9, 249.9, 86, 231.4, 
    42.9, 15.7, 10, 77.8, 273.6, 29.3, 1.5, 86.1, 483.4, 84.8, 
    53.4, 50.3, 131.3, 110.2, 42.5, 9.45, 3.79, 67.8, 97.1, 10.7, 
    144.4, 62.6, 171.6, 29.2, 36.2, 179.6, 383.5, 358.4, 25, 
    4.86, 2.39, 32.4, 75.2, 7.85, 436.1, 630.1, 57.2, 27.7, 25.8, 
    73.7, 314.5, 141.1, 20.2, 1.9, 175.7, 22.8, 68, 4.82, 329.9, 
    573.8, 34.5, 55.1, 40.1, 150.3, 416.7, 44.7, 22.2, 0.87, 
    267.3, 37, 65.3, 5.41, 533.9, 255.2, 69.3, 36.8, 158.8, 95.7, 
    271.1, 26.2, 15.8, 0.81, 45.3, 19.9, 74.3, 2.15, 222.2, 59.5, 
    36, 29.5, 38.7, 56.4, 65.2, 24.8, 29.2, 0.04, 28.4, 17.2, 
    19.5, 1.31, 49, 30.8, 29.6, 39.5, 19.5, 35.8, 123.8, 92.7, 
    31.8, 0, 44.4, 152.2, 17.8, 3.72, 33.4, 44.5, 58.2, 49.4, 
    19.6, 34.4, 692.2)), class = "data.frame", row.names = c(NA, 
-336L))

What I'm trying to do is categorize the different discharge rates into high, medium, and low discharge by month and year.

However, I'm not sure of the most statistically valid way to do this that would make the most sense.

I tried breaking these up into quantiles, but I don't know if that's really the best way to look at it?

I want to use these discharge groups as a factor for analyzing my fish catch data between different months to see how changes in discharge rates effects the fish communities.

$\endgroup$
2
  • 3
    $\begingroup$ What does "river location" mean? Are these different rivers, or are these different spots in a river (eg, just before largest tributary joins, where empties into ocean)? Why not just use the continuous variable as a measure of the background distribution? Do you have reason to believe there are latent groupings (clusters)? Do you just wonder if there are? $\endgroup$ – gung - Reinstate Monica Jul 24 '20 at 19:33
  • 1
    $\begingroup$ Yes, these are two different rivers (B and D) I edited the post to make it more clear. I'm looking to use this as a factor for creating an MDS plot and other analysis in PRIMER. Yes, looking at the data it is clear that we had some drought years and some wetter years that caused a variation in river discharge and I would like to see how these changes impacted the fish communities. Hope this helps to make some sense. $\endgroup$ – fishchick Jul 24 '20 at 19:56
4
$\begingroup$

Here is some basic exploratory data analyses with your data. First, I read in your data and clean it up a little. In particular, I make a new variable, elapsed.months, that is a count of how many months since the beginning of your dataset.

d = structure(list(ï..Zone = c("D", "D", "D", "D", "D", "D", "D",
              ...
              19.6, 34.4, 692.2)), class = "data.frame", row.names = c(NA, -336L))
d$Month          = factor(d$Month, levels=c("Jan", "Feb", "Mar", "April", "May", 
                                            "June", "July", "August", "September", 
                                            "October", "November", "December"))
d$Month          = as.numeric(d$Month)  # this turns them into 1-12
names(d)[1]      = "Zone"
d$Zone           = factor(d$Zone)
d                = d[with(d, order(Zone, Year, Month)),]  # sorted by date
rownames(d)      = NULL
d$elapsed.months = ((d$Year-2005)*12) + d$Month  # months since starting

summary(d$Discharge)  # huge range, minimum is no discharge at all
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#     0.0    25.6    66.9   194.6   249.9  1637.0 
sum(d$Discharge==0)      # [1] 1  # only 1 month without discharge
head(sort(d$Discharge))  # [1] 0.00 0.04 0.13 0.81 0.87 1.31  # lowest 6 values

windows()
  layout(matrix(c(1,1,2,3), nrow=2, byrow=T))
  plot(Discharge~elapsed.months, d, subset=d$Zone=="B", type="l", ylim=c(0,1650))
  with(d[d$Zone=="D",], lines(elapsed.months, Discharge, col="red"))
  legend("topleft", legend=c("D","B"), col=c("red","black"), lty=1)

  boxplot(Discharge~Zone*Month, d, axes=F, col=rep(c("gray","red"), times=12),
          xlab="", main="Months")
  axis(side=2, at=seq(0,1500, by=500));  box()
  axis(side=1, at=seq(0.5, 24.5, by=2), labels=F)
  axis(side=1, at=seq(1.5,23.5, by=2), cex.axis=.8, tick=F, 
       labels=c("J","F","M","A","M","J","J","A","S","O","N","D"))

  boxplot(Discharge~Zone*Year, d, col=rep(c("gray","red"), times=12), axes=F,
          xlab="", main="Months")
  axis(side=2, at=seq(0,1500, by=500));  box()
  axis(side=1, at=seq(0.5, 29.5, by=2), labels=F)
  axis(side=1, at=seq(1.5, 25.5, by=4), cex.axis=.8, tick=F,
       labels=c("05","07","09","11","13","15","17"))
  axis(side=1, at=seq(3.5, 27.5, by=4), cex.axis=.8, tick=F, line=1, 
       labels=c("06","08","10","12","14","16","18"))

line plot of flows over time, and boxplots of flows by month and by year, each stratified by Zone

In the above plots, I see that the two rivers follow a similar pattern over time. River D has higher values and is more variable. That pattern, where variability scales with the mean, is common in data. It often implies a transformation would be appropriate. We can use qq-plots to see if the distributions are roughly normal, and if not, what transformation might be helpful.

library(MASS)  # we'll use this package for the Box-Cox analysis
windows()
  layout(matrix(1:4, nrow=2, byrow=T))
  qqnorm(d$Discharge)
  qqline(d$Discharge)  # nope, not very normal
  bc = with(d, boxcox((Discharge+1)~Month*Year*Zone))
  bc$x[which.max(bc$y)]  # [1] 0.06060606  # suggested optimal transformation
  qqnorm((d$Discharge+1)^.06)
  qqline((d$Discharge+1)^.06)
  qqnorm(log2(d$Discharge+1))
  qqline(log2(d$Discharge+1))

qq-plots and plot of log-likelihoods of Box-Cox parameter

The first qq-plot shows that the data are not remotely normal. Then I use a Box-Cox analysis to find the optimal transformation towards normality. The 95% confidence interval is really narrow and doesn't quite include 0. It suggests the optimal transformation is to take the data to the $.06$ power. However, it is typically recommended to use the nearest interpretable transformation, which in this case is the log (i.e., $0$). I make additional qq-plots for the two transformations, and the difference is negligible to my eye. So we can try using the log of Discharge going forward. Note that I add $1$ to all the original data to avoid taking the log of $0$, and that I use log base $2$ for the transformed data so that each increment of the new values represents a doubling of the flow. The fact that logs work well implies that discharge levels are the result of contributing factors that combine in a multiplicative, rather than additive, manner. Normality isn't quite achieved, though. We have the appearance of bounds on the possible rates. The floor effect is obviously because the rates can't go below 0. The ceiling effect implies there is some constraint that limits the maximum flow possible. Someone who knows more about geology and ecology than I do might be able to advise.

Now we can remake our plots above with the transformed data (code omitted).

d$log.discharge = log2(d$Discharge + 1)
noquote(t(aggregate(log.discharge~Zone, d, summary)))
# Zone                  B         D        
# log.discharge.Min.     0.000000  1.526069
# log.discharge.1st Qu.  4.888729  4.323563
# log.discharge.Median   5.723211  6.666700
# log.discharge.Mean     5.812461  6.409710
# log.discharge.3rd Qu.  7.342219  8.467092
# log.discharge.Max.     9.437128 10.677720
## MADM: the median absolute difference from the median
aggregate(log.discharge~Zone, d, function(x){  median(abs(x-median(x)))  })
#   Zone log.discharge
# 1    B      1.237367
# 2    D      1.975706
1.975706/1.237367  # [1] 1.596702

line and box plots with log transformed data

This looks nicer to me. The values for D continue to be higher than B (the log transformation is monotonic), but the difference is now trivial. On the other hand, D continues to be more variable. Some interesting patterns are discernible. Looking over the years, there was a roughly symmetrical, concave up, parabola between 2005 and 2010. It started downward again, and then abruptly came up to a higher and relatively stable level in 2012-13, where it stayed thereafter. Looking at the months, May through September is much more variable than the rest of the year. October through December typically have lower discharge rates.

Your question is about whether there are meaningful latent groupings. That is ultimately a question for cluster analysis. This case is one-dimensional clustering, but it is still done. A common approach would be to use finite Gaussian mixture modeling, but I wouldn't be comfortable with that here, given the seeming bounded nature of the distribution. Here, I'll just look at kernel density plots of the distribution. We could look at the marginal (overall) distribution, but that doesn't seem like a good idea, because we know there are changes due to seasonality and trends over time. I will use a nonparametric LOWESS smoother to extract the seasonality and trend, and examine the residuals.

ts.D  = ts(d$log.discharge[d$Zone=="D"], frequency=12)
ts.B  = ts(d$log.discharge[d$Zone=="B"], frequency=12)
stl.D = stl(ts.D, "periodic")
stl.B = stl(ts.B, "periodic")
windows()
  plot(stl.D, main="Zone D")
  plot(stl.B, main="Zone B")
windows()
  layout(matrix(1:4, nrow=2, byrow=T))
  plot(density(d$log.discharge[d$Zone=="D"]), main="Marginal distribution for D")
  plot(density(d$log.discharge[d$Zone=="B"]), main="Marginal distribution for B")
  plot(density(stl.D$time.series[,3]),        main="Residual distribution for D")
  plot(density(stl.B$time.series[,3]),        main="Residual distribution for B")

enter image description here
enter image description here
enter image description here

I really don't see anything here to justify binning the data into low, medium and high, or into quartiles, or any other binning scheme. I would just use the discharge rates as a continuous covariate to control for that effect. Note that covariates don't actually have to be normally distributed, contrary to a common belief. Whether you use the raw discharges or the log transformed version is ultimately a theoretical decision based on what makes substantive sense. You could also include both, you only lose one degree of freedom and there's no problem if one isn't significant.

$\endgroup$
3
$\begingroup$

Is there a particular reason that you want to create factors? Doing so will cost you degrees of freedom and loss of information. What is your dependent variable? I would make use of a Month*Discharge interaction to model the discharge effect by month, where the type of regression model will depend on the dependent variable. Make sure to include the first-order versions of Month and Discharge in your model.

*Note - in your sample code, Year is currently coded as numeric. You probably want to convert it to a factor:

df$Year = as.factor(df$Year)
$\endgroup$
2
  • 5
    $\begingroup$ +1 There is extensive discussion here to support your main argument. $\endgroup$ – EdM Jul 24 '20 at 19:59
  • 1
    $\begingroup$ Thank you for your reply and thanks for the note about year. So, I am trying to use the categories of high, medium, low discharge rates as factors for my species catch data to create an MDS plot (among other things) to see if there are trends in fish communities between months or years for each river (zone). I didn't include my catch data because I was trying not to bog the question down too much, but I can add that in if it would help? I apologize for my lack of clarity, I am not great with statistics and I'm still learning these programs. $\endgroup$ – fishchick Jul 24 '20 at 20:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.