Heatmap with points and smoothing for mesocosm experiment. Position of points is (salinity, pH), value/coloring is oogonia.

x=unit_avg_sal y=unit_ph_c x=avg_disc_oo

Helpful sources https://www.r-graph-gallery.com/201-levelplot-with-latticeextra.html
https://latticeextra.r-forge.r-project.org/man/panel.2dsmoother.html
https://www.rdocumentation.org/packages/lattice/versions/0.10-10/topics/panel.levelplot https://www.rdocumentation.org/packages/lattice/versions/0.15-3/topics/levelplot
https://rdrr.io/cran/lattice/man/levelplot.html
https://stackoverflow.com/questions/3712402/r-how-to-change-lattice-levelplot-color-theme #coloring

rm(list=ls())

library(tidyverse)
library(ggpubr)
library(scales)
library(chron)
library(plotly)
library(taRifx)
library(aweek)
library(easypackages)
library(renv)
library(here)
library(ggthemes)
library(gridExtra)
library(patchwork)
library(tidyquant)
library(recipes) 
library(cranlogs)
library(knitr)
library(openair)

#new
library(latticeExtra)

Data

#read in data
heatmap<-read.csv("C:/Users/chels/Box Sync/Thesis/Data/Working data/Final data/meso.calculated.csv", header = TRUE, sep=",", fileEncoding="UTF-8-BOM", stringsAsFactors = FALSE)

#remove rows with NA in acg_disc_ooc
heatmap<-heatmap[!is.na(heatmap$avg_disc_oo), ]

Heatmap

# showing data points on the same color scale 
levelplot(
  avg_disc_oo ~ unit_avg_sal * unit_avg_ph_c,
  heatmap,
  xlab = "Salinity",
  ylab = "pH",
  main = "Effect of salinity and pH on oogonia release",
  panel = panel.levelplot.points,
  col.regions = NULL,
  contour = FALSE,
  cex = 1
) +
  layer_(panel.2dsmoother(..., n = 1000))

####Splitting heatmap to low and high salinity####
Low salinity under 28, high salinity over 30

low.sal<-subset(heatmap, unit_avg_sal <= 28)
hi.sal<-subset(heatmap, unit_avg_sal >= 30)

https://stackoverflow.com/questions/16847377/how-to-set-the-maximum-and-minimum-values-for-a-lattice-plot

Heatmap low salinity

#custom set z range so that the colors of the spilt map can be compared
rgb.palette <- colorRampPalette(c("red", "orange", "yellow", "green", "blue"),
                             space = "rgb")

low.plot<-
  levelplot(
  avg_disc_oo ~ unit_avg_sal * unit_avg_ph_c,
  low.sal,
  xlab = "Low salinity treatments",
  ylab = "pH",
  panel = panel.levelplot.points,
  col.regions = rgb.palette(16),
  at = seq(0, 1800, length=15),
  contour = TRUE,
  colorkey = list(at = seq(0, 1800, length=15),
                  col = rgb.palette(16))) +
    layer_(panel.2dsmoother(..., n = 1000))  

low.plot

Heatmap of high salinity

# showing data points on the same color scale 
hi.plot<-
levelplot(
  avg_disc_oo ~ unit_avg_sal * unit_avg_ph_c,
  hi.sal,
  xlab = "High salinity treatments",
  ylab = "pH",
  panel = panel.levelplot.points,
  col.regions = rgb.palette(16),
  at = seq(0, 1800, length=15),
  contour = TRUE,
  colorkey = list(at = seq(0, 1800, length=15),
                  col = rgb.palette(16))) +
    layer_(panel.2dsmoother(..., n = 1000))

hi.plot

ggarrange

two<-ggarrange(low.plot, hi.plot, ncol=2, nrow=1)

two+plot_annotation(title='Effect of salinity and pH on oogonia release',
                    subtitle='Data from mesocosm experiment',
                    theme = theme(plot.title = element_text(size = 30, family="Calibri Light"),
                                  plot.subtitle = element_text(size=20, family="Calibri Light")))
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database

## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

Need to set the y axis to the same scale
####Looking at linear model####
Heatmap with simple lm: LM not a good fit

levelplot(avg_disc_oo ~ unit_avg_sal * unit_avg_ph_c, heatmap, 
          panel = panel.levelplot.points, cex = 1.2
    ) + 
    layer_(panel.2dsmoother(..., n = 1000,  method = "lm"))

Linear model all data - not a good fit

lm1 <- lm(avg_disc_oo ~ unit_avg_sal + unit_avg_ph_c + unit_avg_sal:unit_avg_ph_c, data=heatmap)  
lm1
## 
## Call:
## lm(formula = avg_disc_oo ~ unit_avg_sal + unit_avg_ph_c + unit_avg_sal:unit_avg_ph_c, 
##     data = heatmap)
## 
## Coefficients:
##                (Intercept)                unit_avg_sal  
##                   3713.960                    -122.592  
##              unit_avg_ph_c  unit_avg_sal:unit_avg_ph_c  
##                   -227.848                       8.293
summary (lm1)
## 
## Call:
## lm(formula = avg_disc_oo ~ unit_avg_sal + unit_avg_ph_c + unit_avg_sal:unit_avg_ph_c, 
##     data = heatmap)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -409.79 -156.97  -73.53  111.96 1264.08 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 3713.960  16500.829   0.225    0.823
## unit_avg_sal                -122.592    602.668  -0.203    0.840
## unit_avg_ph_c               -227.848   2098.387  -0.109    0.914
## unit_avg_sal:unit_avg_ph_c     8.293     76.529   0.108    0.914
## 
## Residual standard error: 319.8 on 56 degrees of freedom
## Multiple R-squared:  0.1299, Adjusted R-squared:  0.08329 
## F-statistic: 2.787 on 3 and 56 DF,  p-value: 0.04897
anova (lm1)
## Analysis of Variance Table
## 
## Response: avg_disc_oo
##                            Df  Sum Sq Mean Sq F value   Pr(>F)   
## unit_avg_sal                1  854078  854078  8.3491 0.005479 **
## unit_avg_ph_c               1       6       6  0.0001 0.994114   
## unit_avg_sal:unit_avg_ph_c  1    1201    1201  0.0117 0.914094   
## Residuals                  56 5728545  102295                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot (lm1)

Linear model looking at high and low salinity separately Low sal- looks like a much better fit

lm2 <- lm(avg_disc_oo ~ unit_avg_sal + unit_avg_ph_c + unit_avg_sal:unit_avg_ph_c, data=low.sal)  
lm2
## 
## Call:
## lm(formula = avg_disc_oo ~ unit_avg_sal + unit_avg_ph_c + unit_avg_sal:unit_avg_ph_c, 
##     data = low.sal)
## 
## Coefficients:
##                (Intercept)                unit_avg_sal  
##                   111520.2                     -4117.4  
##              unit_avg_ph_c  unit_avg_sal:unit_avg_ph_c  
##                   -13601.4                       503.5
summary (lm2)
## 
## Call:
## lm(formula = avg_disc_oo ~ unit_avg_sal + unit_avg_ph_c + unit_avg_sal:unit_avg_ph_c, 
##     data = low.sal)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -552.83 -266.18  -87.47  193.84 1093.21 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                111520.2   116237.0   0.959    0.346
## unit_avg_sal                -4117.4     4317.8  -0.954    0.349
## unit_avg_ph_c              -13601.4    15285.8  -0.890    0.382
## unit_avg_sal:unit_avg_ph_c    503.5      568.1   0.886    0.384
## 
## Residual standard error: 414.9 on 26 degrees of freedom
## Multiple R-squared:  0.07898,    Adjusted R-squared:  -0.02729 
## F-statistic: 0.7432 on 3 and 26 DF,  p-value: 0.536
anova (lm2)
## Analysis of Variance Table
## 
## Response: avg_disc_oo
##                            Df  Sum Sq Mean Sq F value Pr(>F)
## unit_avg_sal                1  223243  223243  1.2968 0.2652
## unit_avg_ph_c               1   25351   25351  0.1473 0.7043
## unit_avg_sal:unit_avg_ph_c  1  135236  135236  0.7856 0.3836
## Residuals                  26 4475993  172154
plot (lm2)

High sal- looks like a better fit but not sure if its good enough

lm3 <- lm(avg_disc_oo ~ unit_avg_sal + unit_avg_ph_c + unit_avg_sal:unit_avg_ph_c, data=hi.sal)  
lm3
## 
## Call:
## lm(formula = avg_disc_oo ~ unit_avg_sal + unit_avg_ph_c + unit_avg_sal:unit_avg_ph_c, 
##     data = hi.sal)
## 
## Coefficients:
##                (Intercept)                unit_avg_sal  
##                  -232053.8                      7610.5  
##              unit_avg_ph_c  unit_avg_sal:unit_avg_ph_c  
##                    27662.2                      -906.9
summary (lm3)
## 
## Call:
## lm(formula = avg_disc_oo ~ unit_avg_sal + unit_avg_ph_c + unit_avg_sal:unit_avg_ph_c, 
##     data = hi.sal)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -204.35  -96.33  -32.70   54.45  516.26 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)                -232053.8   363972.8  -0.638    0.529
## unit_avg_sal                  7610.5    11883.8   0.640    0.528
## unit_avg_ph_c                27662.2    45937.3   0.602    0.552
## unit_avg_sal:unit_avg_ph_c    -906.9     1499.8  -0.605    0.551
## 
## Residual standard error: 176.4 on 26 degrees of freedom
## Multiple R-squared:  0.1302, Adjusted R-squared:  0.02984 
## F-statistic: 1.297 on 3 and 26 DF,  p-value: 0.2964
anova (lm3)
## Analysis of Variance Table
## 
## Response: avg_disc_oo
##                            Df Sum Sq Mean Sq F value  Pr(>F)  
## unit_avg_sal                1  91537   91537  2.9430 0.09815 .
## unit_avg_ph_c               1  18144   18144  0.5833 0.45188  
## unit_avg_sal:unit_avg_ph_c  1  11372   11372  0.3656 0.55063  
## Residuals                  26 808687   31103                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot (lm3)