search

andrewcparnell / simmr

A stable isotope mixing model in R
https://andrewcparnell.github.io/simmr/
28 stars 8 forks source link

Changing limits on Y axis for compare_groups #28

Closed Craigdux closed 3 years ago

Craigdux commented 3 years ago

Hi Andrew, Your mixing model simmr is a huge help in our work

I have outputted some nice box plots with the compare_groups function. However, I cant figure out how to change the default title, nor the limit on the y-axis.

I am familiar with ggplot routines, but I can't seem to feed the plot output into ggplot.

Can you offer any advice?
Thanks! Craig Duxbury

andrewcparnell commented 3 years ago

Hi Craig,

Thanks for using simmr. Have you seen this section of the manual? Let me know if it helps at all.

Andrew

Craigdux commented 3 years ago

Hi Andrew--Thanks for your reply. My problem was (I am fairly new user of R) is that I was not sure how to call up the plot from the compare_groups function.

I had someone help me by first: Plot1 <- compare_groups....

Then Plot1$plot + labs(x = "Well", title = "Denitrification).

Worked like a charm!

Again thanks, and it was my lack of understanding of the output from your compare_groups.

Best,

Craig

andrewcparnell commented 3 years ago

Well done! Let me know if there's anything that needs improving.

Andrew

Craigdux commented 3 years ago

I do have one more question:

I was able to make a convex hull with simmr, however, I was not able to superimpose source areas, similar to what I have here. I am inheriting this analyses from someone who did it in siar, hence my problem: image

andrewcparnell commented 3 years ago

Not sure exactly what you mean. Are you saying you couldn't produce that plot (done in PowerPoint or similar)? Or that you were trying to produce something with different areas on top of it?

Craigdux commented 3 years ago

I created the convex hull plot in simmr and the source plot separately in simmr. But, as you noted, I had to use Powerpoint to draw the convex hull on the source plot.

Is there a way to combine them?

thank you.

andrewcparnell commented 3 years ago

OK I've written a new vignette on how to customise plots which includes a way of adding convex hulls and polygons on to plots. See here: https://andrewcparnell.github.io/simmr/articles/advanced_plotting.html

Hope it helps!

Andrew

Craigdux commented 3 years ago

@andrewcparnell --Thank you for doing this!!!

Craigdux commented 2 years ago

@andrewcparnell --sorry to be a pain in your side on this one, but I am trying to get this plot completed and it is throwing an error when I try to add the convex hull. I copied your code and created a biplot. Then using your code from the vignette, and using my original "sources" data frame for the convex hull. When I run the "chull" command, it does not change the values for plotting (that does make sense), It seems I am doing something simple--perhaps I should not be using my "sources". Can you help? Thank you.

The error:

Error in check_aesthetics(): ! Aesthetics must be either length 1 or the same as the data (4): colour_new

My code:

Isoplot Biplot

p<- plot(well_simmr, group = 1:8, xlab = expression(paste(delta^15, "N (\u2030)", sep = "")), ylab = expression(paste(delta^18, "O (\u2030)", sep = "")), title = 'Test', mix_name = '')

p + ggtitle("Test2")

Add convex hulls

chull_vals <- chull(sources[ ,2], sources[ , 3]) sources2 <- sources[chull_vals, ]

p + new_scale_color() + geom_polygon(data = sources2, aes(x = meand15N, y = meand18O), fill = "orange", alpha = 0.2)

And the biplot: image

andrewcparnell commented 2 years ago

I don't think I can tell from that code. Can you make a reproducible example?

Craigdux commented 2 years ago

OK. I will add my pertinent code (the simmr model parts are running amazingly, BTW!!--Thank you). It is just this polygon that is giving me grief! Please let me know if you can run this. Thank you

The data

data.frame(
  stringsAsFactors = FALSE,
                                   d15N = c(4.6,2,2.9,6.3,4.1,
                                            2.9,4.4,9.7,4.3,2.3,3.3,8.4,5.5,
                                            2.6,4.5,4.7,6.9,4.5,7.7,2.8,
                                            6.3,4.2,7.9,6.4,3.9,3.6,2.5,
                                            1.9,6.4,2.1,2.7,2.5,2.8,2.5,
                                            2.5,3,3.3,2.9,6.5,2.4,2.8,6.2,
                                            7.7,2.9,13.4,5.8,3,3.1,9.7,
                                            4.4,5.4,3.8,5.7,4.8,2.9,2.7,
                                            6.8,6.4,3.3,6.4,2.1,8.3,4.6,3.1,
                                            6.2,6.7,2.8,8.9,3.9,2.6,5.7,
                                            5.9,7.1,5.1,9.6,2.9,3,7.3,
                                            8.3,7.2,5.8,9.1,10.7,5.1,12.8,
                                            5.4,7.4,6.1,2.6,7.4,10.9,4.7,
                                            6.9,3,5.1,7.1,3.8,10.7,6.3,3,
                                            5.9,6.7,7.1,5.4,2.4,9.5,10.7,
                                            5.1,4.7,3,6.7,6.9,1.9,8.2,
                                            11.7,4.5,5.8,2.9,7.3,6.4,1.7,7,
                                            13,8.7,8.1,2.7,8.5,7.5,3.5,9,
                                            14.5,8.4,11.5,3.7,6.4,7.7,3.4,
                                            10.2,7.2,4.3,12.6,8.3,10.9,
                                            6.2,5.8,10.8,10,6.7,8,6.1,10.2,
                                            6.3,5.8,10.4,8.3,6.4,8.2,5.6,
                                            8.3,7.5,6,7.7,9.5,5.9,5.8,
                                            5.5),
                                   d180 = c(3.5,2,2.3,4.7,3.2,
                                            1.2,4,11.4,10.5,4.9,2.1,6,4.2,
                                            1.7,2.6,1,5.5,2.1,10.2,2.7,4,
                                            2.5,4.9,7.4,5.6,4.4,2.4,1.8,
                                            11.3,11.3,3.4,3.4,2.4,7.8,2.2,
                                            4.4,4.8,1.7,11,-0.1,2.2,9.7,
                                            9.9,-2.5,16,3.9,1.8,6,12.2,
                                            3.3,4.8,1.1,4.5,4.8,8,4.8,6.8,
                                            12.6,2.2,12.2,5.3,8.9,7.8,3.3,
                                            2.8,15.4,12.4,15.3,2.1,2.5,
                                            4.3,7.1,12.7,7.6,8.1,-1.9,0.1,
                                            8.1,14.1,13.3,3.5,13.6,9.2,7,
                                            18.8,2.9,12.4,10.1,7,17.3,15.8,
                                            6.9,9.6,5.1,7,12.1,2.9,10.7,
                                            6.6,2.6,7.1,8.8,4.4,3,1.1,5.6,
                                            14.3,13.5,7.6,3.6,6.4,4.6,
                                            1.3,12.8,11.3,5,5.2,-0.1,7.5,
                                            4.5,2,12,19.6,22.7,6.7,1.2,6.6,
                                            3.2,1.7,11.5,20.8,9.7,7.9,2.3,
                                            20.9,12.6,9.1,20.2,6.4,2.9,
                                            16.6,15.9,7.2,5.7,3.4,12.3,7,
                                            5.9,8.1,2.8,7.7,4.3,3.6,13,
                                            10.3,3.7,7.8,2.4,6.1,4.9,3.3,
                                            10.8,12.2,4.7,6.7,3.2),
                                   well = c("MW22","MW07","MW04",
                                            "MW11","MW22","MW07","MW11",
                                            "MW20","MW22","MW07","MW11","MWBS",
                                            "MW22","MW07","MW11","MW07",
                                            "MW11","MW04","MWDS","MW07","MW22",
                                            "MW04","MW11","MW22","MW07",
                                            "MW11","MW04","MW22","MWDS","MW07",
                                            "MW11","MW04","MW11","MW22",
                                            "MW04","MW11","MW22","MW04","MWDS",
                                            "MW07","MW04","MWDS","MWDS",
                                            "MW07","MW17","MW11","MW04","MW22",
                                            "MWDS","MW07","MW11","MW04",
                                            "MWBS","MW22","MW07","MW11","MWBS",
                                            "MWDS","MW04","MWDS","MW22",
                                            "MW17","MW07","MW04","MWBS","MWDS",
                                            "MW22","MW17","MW07","MW04",
                                            "MW11","MWBS","MWDS","MW22","MW17",
                                            "MW07","MW04","MWBS","MWBS",
                                            "MW11","MW04","MWDS","MW17","MW22",
                                            "MW20","MW07","MWBS","MW11",
                                            "MW04","MWDS","MW17","MW22","MW20",
                                            "MW07","MW22","MWDS","MW07",
                                            "MW17","MW11","MW04","MW20","MWBS",
                                            "MWBS","MW11","MW04","MWDS",
                                            "MW17","MW22","MW20","MW07","MWBS",
                                            "MW11","MW04","MWDS","MW17",
                                            "MW22","MW20","MW07","MWBS","MW11",
                                            "MW04","MWDS","MW17","MW22",
                                            "MW20","MW07","MWBS","MW11","MW04",
                                            "MWDS","MW17","MW22","MW20",
                                            "MW07","MW22","MWDS","MW07","MW17",
                                            "MW11","MW04","MW20","MWBS",
                                            "MWBS","MW11","MW04","MWDS","MW17",
                                            "MW22","MW20","MW07","MWBS",
                                            "MW11","MW04","MWDS","MW17","MW22",
                                            "MW20","MW07","MWBS","MW11",
                                            "MW04","MWDS","MW17","MW22","MW20",
                                            "MW07")
                     )
#>     d15N d180 well
#> 1    4.6  3.5 MW22
#> 2    2.0  2.0 MW07
#> 3    2.9  2.3 MW04
#> 4    6.3  4.7 MW11
#> 5    4.1  3.2 MW22
#> 6    2.9  1.2 MW07
#> 7    4.4  4.0 MW11
#> 8    9.7 11.4 MW20
#> 9    4.3 10.5 MW22
#> 10   2.3  4.9 MW07
#> 11   3.3  2.1 MW11
#> 12   8.4  6.0 MWBS
#> 13   5.5  4.2 MW22
#> 14   2.6  1.7 MW07
#> 15   4.5  2.6 MW11
#> 16   4.7  1.0 MW07
#> 17   6.9  5.5 MW11
#> 18   4.5  2.1 MW04
#> 19   7.7 10.2 MWDS
#> 20   2.8  2.7 MW07
#> 21   6.3  4.0 MW22
#> 22   4.2  2.5 MW04
#> 23   7.9  4.9 MW11
#> 24   6.4  7.4 MW22
#> 25   3.9  5.6 MW07
#> 26   3.6  4.4 MW11
#> 27   2.5  2.4 MW04
#> 28   1.9  1.8 MW22
#> 29   6.4 11.3 MWDS
#> 30   2.1 11.3 MW07
#> 31   2.7  3.4 MW11
#> 32   2.5  3.4 MW04
#> 33   2.8  2.4 MW11
#> 34   2.5  7.8 MW22
#> 35   2.5  2.2 MW04
#> 36   3.0  4.4 MW11
#> 37   3.3  4.8 MW22
#> 38   2.9  1.7 MW04
#> 39   6.5 11.0 MWDS
#> 40   2.4 -0.1 MW07
#> 41   2.8  2.2 MW04
#> 42   6.2  9.7 MWDS
#> 43   7.7  9.9 MWDS
#> 44   2.9 -2.5 MW07
#> 45  13.4 16.0 MW17
#> 46   5.8  3.9 MW11
#> 47   3.0  1.8 MW04
#> 48   3.1  6.0 MW22
#> 49   9.7 12.2 MWDS
#> 50   4.4  3.3 MW07
#> 51   5.4  4.8 MW11
#> 52   3.8  1.1 MW04
#> 53   5.7  4.5 MWBS
#> 54   4.8  4.8 MW22
#> 55   2.9  8.0 MW07
#> 56   2.7  4.8 MW11
#> 57   6.8  6.8 MWBS
#> 58   6.4 12.6 MWDS
#> 59   3.3  2.2 MW04
#> 60   6.4 12.2 MWDS
#> 61   2.1  5.3 MW22
#> 62   8.3  8.9 MW17
#> 63   4.6  7.8 MW07
#> 64   3.1  3.3 MW04
#> 65   6.2  2.8 MWBS
#> 66   6.7 15.4 MWDS
#> 67   2.8 12.4 MW22
#> 68   8.9 15.3 MW17
#> 69   3.9  2.1 MW07
#> 70   2.6  2.5 MW04
#> 71   5.7  4.3 MW11
#> 72   5.9  7.1 MWBS
#> 73   7.1 12.7 MWDS
#> 74   5.1  7.6 MW22
#> 75   9.6  8.1 MW17
#> 76   2.9 -1.9 MW07
#> 77   3.0  0.1 MW04
#> 78   7.3  8.1 MWBS
#> 79   8.3 14.1 MWBS
#> 80   7.2 13.3 MW11
#> 81   5.8  3.5 MW04
#> 82   9.1 13.6 MWDS
#> 83  10.7  9.2 MW17
#> 84   5.1  7.0 MW22
#> 85  12.8 18.8 MW20
#> 86   5.4  2.9 MW07
#> 87   7.4 12.4 MWBS
#> 88   6.1 10.1 MW11
#> 89   2.6  7.0 MW04
#> 90   7.4 17.3 MWDS
#> 91  10.9 15.8 MW17
#> 92   4.7  6.9 MW22
#> 93   6.9  9.6 MW20
#> 94   3.0  5.1 MW07
#> 95   5.1  7.0 MW22
#> 96   7.1 12.1 MWDS
#> 97   3.8  2.9 MW07
#> 98  10.7 10.7 MW17
#> 99   6.3  6.6 MW11
#> 100  3.0  2.6 MW04
#> 101  5.9  7.1 MW20
#> 102  6.7  8.8 MWBS
#> 103  7.1  4.4 MWBS
#> 104  5.4  3.0 MW11
#> 105  2.4  1.1 MW04
#> 106  9.5  5.6 MWDS
#> 107 10.7 14.3 MW17
#> 108  5.1 13.5 MW22
#> 109  4.7  7.6 MW20
#> 110  3.0  3.6 MW07
#> 111  6.7  6.4 MWBS
#> 112  6.9  4.6 MW11
#> 113  1.9  1.3 MW04
#> 114  8.2 12.8 MWDS
#> 115 11.7 11.3 MW17
#> 116  4.5  5.0 MW22
#> 117  5.8  5.2 MW20
#> 118  2.9 -0.1 MW07
#> 119  7.3  7.5 MWBS
#> 120  6.4  4.5 MW11
#> 121  1.7  2.0 MW04
#> 122  7.0 12.0 MWDS
#> 123 13.0 19.6 MW17
#> 124  8.7 22.7 MW22
#> 125  8.1  6.7 MW20
#> 126  2.7  1.2 MW07
#> 127  8.5  6.6 MWBS
#> 128  7.5  3.2 MW11
#> 129  3.5  1.7 MW04
#> 130  9.0 11.5 MWDS
#> 131 14.5 20.8 MW17
#> 132  8.4  9.7 MW22
#> 133 11.5  7.9 MW20
#> 134  3.7  2.3 MW07
#> 135  6.4 20.9 MW22
#> 136  7.7 12.6 MWDS
#> 137  3.4  9.1 MW07
#> 138 10.2 20.2 MW17
#> 139  7.2  6.4 MW11
#> 140  4.3  2.9 MW04
#> 141 12.6 16.6 MW20
#> 142  8.3 15.9 MWBS
#> 143 10.9  7.2 MWBS
#> 144  6.2  5.7 MW11
#> 145  5.8  3.4 MW04
#> 146 10.8 12.3 MWDS
#> 147 10.0  7.0 MW17
#> 148  6.7  5.9 MW22
#> 149  8.0  8.1 MW20
#> 150  6.1  2.8 MW07
#> 151 10.2  7.7 MWBS
#> 152  6.3  4.3 MW11
#> 153  5.8  3.6 MW04
#> 154 10.4 13.0 MWDS
#> 155  8.3 10.3 MW17
#> 156  6.4  3.7 MW22
#> 157  8.2  7.8 MW20
#> 158  5.6  2.4 MW07
#> 159  8.3  6.1 MWBS
#> 160  7.5  4.9 MW11
#> 161  6.0  3.3 MW04
#> 162  7.7 10.8 MWDS
#> 163  9.5 12.2 MW17
#> 164  5.9  4.7 MW22
#> 165  5.8  6.7 MW20
#> 166  5.5  3.2 MW07

Created on 2022-10-17 by the reprex package (v2.0.1)

The sources:

data.frame(
  stringsAsFactors = FALSE,
                      Sources = c("NH4 Fertilizer","Septic and Manure","NO3 Fertilizer",
                                  "Denitrification"),
          meand15N = c(-0.9, 12.8, 0.65, 25),
          meand18O = c(4.18, 4.18, 8.54, 25),
            SDd15N = c(2.07, 6.1, 1.75, 5),
            SDd18O = c(2.87, 2.87, 2.92, 5)
           )
#>             Sources meand15N meand18O SDd15N SDd18O
#> 1    NH4 Fertilizer    -0.90     4.18   2.07   2.87
#> 2 Septic and Manure    12.80     4.18   6.10   2.87
#> 3    NO3 Fertilizer     0.65     8.54   1.75   2.92
#> 4   Denitrification    25.00    25.00   5.00   5.00

Created on 2022-10-17 by the reprex package (v2.0.1)

My Code

## Bring in data from  csv files
targets <-read.csv("well.isotopes.csv")   ### Actual raw data
#> Warning in file(file, "rt"): cannot open file 'well.isotopes.csv': No such file
#> or directory
#> Error in file(file, "rt"): cannot open the connection
sources <- read.csv("sources.csv")   ### Means and SD of sources
#> Warning in file(file, "rt"): cannot open file 'sources.csv': No such file or
#> directory
#> Error in file(file, "rt"): cannot open the connection
tefs <- read.csv("tefs.csv")
#> Warning in file(file, "rt"): cannot open file 'tefs.csv': No such file or
#> directory
#> Error in file(file, "rt"): cannot open the connection
concep <- read.csv("concep.csv")
#> Warning in file(file, "rt"): cannot open file 'concep.csv': No such file or
#> directory
#> Error in file(file, "rt"): cannot open the connection

### Setup the model
well_simmr <- simmr_load(mixtures = as.matrix(targets[, 1:2]),
                        source_names = sources$Sources,
                        source_means = as.matrix(sources[,2:3]),
                        source_sds = as.matrix(sources[,4:5]),
                        #correction_means = as.matrix(tefs[,2:3]),
                        #correction_sds = as.matrix(tefs[,4:5]),
                        #concentration_means = as.matrix(concep[,2:3]),
                        group = targets$well)
#> Error in simmr_load(mixtures = as.matrix(targets[, 1:2]), source_names = sources$Sources, : could not find function "simmr_load"

### Isoplot Biplot
p<- plot(well_simmr,
         group = 1:8,
         xlab = expression(paste(delta^15, "N (\u2030)",
                                 sep = "")), 
         ylab = expression(paste(delta^18, "O (\u2030)",
                                 sep = "")), 
         title = 'Test',
         mix_name = '')
#> Error in plot(well_simmr, group = 1:8, xlab = expression(paste(delta^15, : object 'well_simmr' not found

## Add convex hulls
chull_vals <- chull(sources[ ,2], sources[ , 3])
#> Error in xy.coords(x, y, recycle = TRUE, setLab = FALSE): object 'sources' not found
sources2 <- sources[chull_vals, ]
#> Error in eval(expr, envir, enclos): object 'sources' not found

p + new_scale_color() +
  geom_polygon(data = sources2, aes(x = meand15N, y = meand18O), fill = "orange", alpha = 0.2)
#> Error in eval(expr, envir, enclos): object 'p' not found

Created on 2022-10-17 by the reprex package (v2.0.1)

andrewcparnell commented 2 years ago

You were missing a few lines in there I think before creating the convex hull. Try this:

# Script for craigdux using simmr
# Having problems with adding a convex hull
library(simmr)
library(ggnewscale)

targets <- data.frame(
  stringsAsFactors = FALSE,
  d15N = c(4.6,2,2.9,6.3,4.1,
           2.9,4.4,9.7,4.3,2.3,3.3,8.4,5.5,
           2.6,4.5,4.7,6.9,4.5,7.7,2.8,
           6.3,4.2,7.9,6.4,3.9,3.6,2.5,
           1.9,6.4,2.1,2.7,2.5,2.8,2.5,
           2.5,3,3.3,2.9,6.5,2.4,2.8,6.2,
           7.7,2.9,13.4,5.8,3,3.1,9.7,
           4.4,5.4,3.8,5.7,4.8,2.9,2.7,
           6.8,6.4,3.3,6.4,2.1,8.3,4.6,3.1,
           6.2,6.7,2.8,8.9,3.9,2.6,5.7,
           5.9,7.1,5.1,9.6,2.9,3,7.3,
           8.3,7.2,5.8,9.1,10.7,5.1,12.8,
           5.4,7.4,6.1,2.6,7.4,10.9,4.7,
           6.9,3,5.1,7.1,3.8,10.7,6.3,3,
           5.9,6.7,7.1,5.4,2.4,9.5,10.7,
           5.1,4.7,3,6.7,6.9,1.9,8.2,
           11.7,4.5,5.8,2.9,7.3,6.4,1.7,7,
           13,8.7,8.1,2.7,8.5,7.5,3.5,9,
           14.5,8.4,11.5,3.7,6.4,7.7,3.4,
           10.2,7.2,4.3,12.6,8.3,10.9,
           6.2,5.8,10.8,10,6.7,8,6.1,10.2,
           6.3,5.8,10.4,8.3,6.4,8.2,5.6,
           8.3,7.5,6,7.7,9.5,5.9,5.8,
           5.5),
  d180 = c(3.5,2,2.3,4.7,3.2,
           1.2,4,11.4,10.5,4.9,2.1,6,4.2,
           1.7,2.6,1,5.5,2.1,10.2,2.7,4,
           2.5,4.9,7.4,5.6,4.4,2.4,1.8,
           11.3,11.3,3.4,3.4,2.4,7.8,2.2,
           4.4,4.8,1.7,11,-0.1,2.2,9.7,
           9.9,-2.5,16,3.9,1.8,6,12.2,
           3.3,4.8,1.1,4.5,4.8,8,4.8,6.8,
           12.6,2.2,12.2,5.3,8.9,7.8,3.3,
           2.8,15.4,12.4,15.3,2.1,2.5,
           4.3,7.1,12.7,7.6,8.1,-1.9,0.1,
           8.1,14.1,13.3,3.5,13.6,9.2,7,
           18.8,2.9,12.4,10.1,7,17.3,15.8,
           6.9,9.6,5.1,7,12.1,2.9,10.7,
           6.6,2.6,7.1,8.8,4.4,3,1.1,5.6,
           14.3,13.5,7.6,3.6,6.4,4.6,
           1.3,12.8,11.3,5,5.2,-0.1,7.5,
           4.5,2,12,19.6,22.7,6.7,1.2,6.6,
           3.2,1.7,11.5,20.8,9.7,7.9,2.3,
           20.9,12.6,9.1,20.2,6.4,2.9,
           16.6,15.9,7.2,5.7,3.4,12.3,7,
           5.9,8.1,2.8,7.7,4.3,3.6,13,
           10.3,3.7,7.8,2.4,6.1,4.9,3.3,
           10.8,12.2,4.7,6.7,3.2),
  well = c("MW22","MW07","MW04",
           "MW11","MW22","MW07","MW11",
           "MW20","MW22","MW07","MW11","MWBS",
           "MW22","MW07","MW11","MW07",
           "MW11","MW04","MWDS","MW07","MW22",
           "MW04","MW11","MW22","MW07",
           "MW11","MW04","MW22","MWDS","MW07",
           "MW11","MW04","MW11","MW22",
           "MW04","MW11","MW22","MW04","MWDS",
           "MW07","MW04","MWDS","MWDS",
           "MW07","MW17","MW11","MW04","MW22",
           "MWDS","MW07","MW11","MW04",
           "MWBS","MW22","MW07","MW11","MWBS",
           "MWDS","MW04","MWDS","MW22",
           "MW17","MW07","MW04","MWBS","MWDS",
           "MW22","MW17","MW07","MW04",
           "MW11","MWBS","MWDS","MW22","MW17",
           "MW07","MW04","MWBS","MWBS",
           "MW11","MW04","MWDS","MW17","MW22",
           "MW20","MW07","MWBS","MW11",
           "MW04","MWDS","MW17","MW22","MW20",
           "MW07","MW22","MWDS","MW07",
           "MW17","MW11","MW04","MW20","MWBS",
           "MWBS","MW11","MW04","MWDS",
           "MW17","MW22","MW20","MW07","MWBS",
           "MW11","MW04","MWDS","MW17",
           "MW22","MW20","MW07","MWBS","MW11",
           "MW04","MWDS","MW17","MW22",
           "MW20","MW07","MWBS","MW11","MW04",
           "MWDS","MW17","MW22","MW20",
           "MW07","MW22","MWDS","MW07","MW17",
           "MW11","MW04","MW20","MWBS",
           "MWBS","MW11","MW04","MWDS","MW17",
           "MW22","MW20","MW07","MWBS",
           "MW11","MW04","MWDS","MW17","MW22",
           "MW20","MW07","MWBS","MW11",
           "MW04","MWDS","MW17","MW22","MW20",
           "MW07")
)

sources <- data.frame(
  stringsAsFactors = FALSE,
  Sources = c("NH4 Fertilizer","Septic and Manure","NO3 Fertilizer",
              "Denitrification"),
  meand15N = c(-0.9, 12.8, 0.65, 25),
  meand18O = c(4.18, 4.18, 8.54, 25),
  SDd15N = c(2.07, 6.1, 1.75, 5),
  SDd18O = c(2.87, 2.87, 2.92, 5)
)

### Setup the model
well_simmr <- simmr_load(mixtures = as.matrix(targets[, 1:2]),
                         source_names = sources$Sources,
                         source_means = as.matrix(sources[,2:3]),
                         source_sds = as.matrix(sources[,4:5]),
                         group = targets$well)

### Isoplot Biplot
p <- plot(well_simmr,
         group = 1:8,
         xlab = expression(paste(delta^15, "N (\u2030)",
                                 sep = "")), 
         ylab = expression(paste(delta^18, "O (\u2030)",
                                 sep = "")), 
         title = 'Test',
         mix_name = '')

## Add convex hulls
new_df <- data.frame(
  x = sources$meand15N,
  y = sources$meand18O,
  Source = levels(p$data$Source)[1] # This just needs to be one of the names of the groups or the sources
)
chull_vals <- chull(new_df[, 1], new_df[, 2])
new_df2 <- new_df[chull_vals, ]
p + new_scale_color() +
  geom_polygon(data = new_df2, aes(x = x, y = y), fill = "orange", alpha = 0.2)

# You can stop here, or...
# Here's a fancier bit if you want it
new_df3 <- data.frame(
  x = c(sources$meand15N - sources$SDd15N, sources$meand15N - sources$SDd15N,
        sources$meand15N + sources$SDd15N, sources$meand15N + sources$SDd15N),
  y = c(sources$meand18O - sources$SDd18O, sources$meand18O + sources$SDd18O,
        sources$meand18O - sources$SDd18O, sources$meand18O + sources$SDd18O),
  Source = levels(p$data$Source)[1] # This just needs to be one of the names of the groups or the sources
)
chull_vals <- chull(new_df3[, 1], new_df3[, 2])
new_df4 <- new_df3[chull_vals, ]
p + new_scale_color() +
  geom_polygon(data = new_df2, aes(x = x, y = y), fill = "orange", alpha = 0.3) + 
  new_scale_color() +
  geom_polygon(data = new_df4, aes(x = x, y = y), fill = "orange", alpha = 0.1)
Craigdux commented 2 years ago

@andrewcparnell This is fantastic! The fancier bit is super helpful, too.

How do interpret the larger convex polygons?

thanks!

andrewcparnell commented 2 years ago

Well if you have mixture points outside the wider polygon it's probably a sign that you're missing a source or they are screwy measurements. Though really these plots only use 1 standard deviation; if you want to produce something more analogous to a 95% source convex hull then they should all be changed to mean +/- 2 standard deviations.

Craigdux commented 2 years ago

Thank you again!!! I owe you a good bottle of wine!
C