library(readxl)
library(lubridate)
library(tidyverse)
library(plotly)
library(stargazer)
library(Metrics)

AIM : To compare HOBO (lux) and LiCOR (PAR) data. Found out conversion rate between LUX and PAR. Explore potential biaises of relationship with Depth or HOBO sensors. Compare linear and non linear regression following Long et al. (2012).

Document Note : The maps, plots and app within this document are interactive so make sure you give them a play like zooming in and out in the maps but also on the plots. Clicking on the legend allows to only select and display the time series needed.

Table of contents

  1. Linear regression PAR vs LUX

  1. Non-linear regression PAR vs LUX

  2. Bibliography

Linear regression PAR vs LUX

All Depth and Sensors blend together 

hobo_licor <- read_excel('Kaik_lux_light_2021.xlsx')

## Plot Hobo vs PAR - No distinction of Depth, Time (hour), Site
p <- ggplot(hobo_licor,aes(x=Lux,y=PAR)) +
  geom_point() + 
  geom_smooth(formula = y ~ x, se=T, method = 'lm') + 
  theme_light() +
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=20))
ggplotly(p)

Figure 1 - PAR vs LUX and linear regression.

\(R^2\) of linear fit:

lm_fit <- lm(PAR ~ Lux, data = hobo_licor) #%>% tidy() %>% kable()
#summary(lm_fit)#$r.squared

#lm(Kd_MW ~ Kd_NN, data = nn_mw_idw_tb) %>% kable()

stargazer(lm_fit,title="Regression Results",type = "text", style = 'default')#"html"
## 
## Regression Results
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                 PAR            
## -----------------------------------------------
## Lux                          0.028***          
##                               (0.001)          
##                                                
## Constant                     68.408***         
##                               (3.031)          
##                                                
## -----------------------------------------------
## Observations                    516            
## R2                             0.817           
## Adjusted R2                    0.817           
## Residual Std. Error      34.468 (df = 514)     
## F Statistic         2,297.848*** (df = 1; 514) 
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

\(R^2\) = .817, probably not too bad but definitely some variability. Due to what? Sensors? Depth?

\(PAR = 0.028*LUX + 68.408\)

Effect of Hobo sensor 

p <- ggplot(hobo_licor,aes(x=Lux,y=PAR,col=Hobo,group=Hobo)) +
  geom_point() + 
  geom_smooth(formula = y ~ x, se=T, method = 'lm') + 
  theme_light() +
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=20))
ggplotly(p)

Figure 2 - PAR vs LUX and effect of sensor type on linear regression.

\(R^2\) of linear fit:

lm_fit_sensor40 <- hobo_licor %>% 
  dplyr::filter(Hobo == '40') %>% 
  lm(PAR ~ Lux, data = .) 

lm_fit_sensorwbn2 <- hobo_licor %>% 
  dplyr::filter(Hobo == 'WBN 2') %>% 
  lm(PAR ~ Lux, data = .) 

lm_fit_sensor5 <- hobo_licor %>% 
  dplyr::filter(Hobo == '5') %>% 
  lm(PAR ~ Lux, data = .) 

lm_fit_sensor21 <- hobo_licor %>% 
  dplyr::filter(Hobo == '21') %>% 
  lm(PAR ~ Lux, data = .) 

stargazer(lm_fit_sensor21,lm_fit_sensor40,lm_fit_sensor5,lm_fit_sensorwbn2, title="Regression Results for Sensors 21 (1), 40 (2), 5 (3) and WBN 2 (4).",type = "text", style = 'default')#"html"
## 
## Regression Results for Sensors 21 (1), 40 (2), 5 (3) and WBN 2 (4).
## ============================================================================
##                                             Dependent variable:             
##                                ---------------------------------------------
##                                                     PAR                     
##                                   (1)        (2)         (3)         (4)    
## ----------------------------------------------------------------------------
## Lux                             0.028***   0.025***    0.033***    0.031*** 
##                                 (0.001)    (0.001)     (0.001)     (0.001)  
##                                                                             
## Constant                       62.953***  66.663***   60.628***   64.244*** 
##                                 (5.160)    (6.538)     (4.723)     (5.813)  
##                                                                             
## ----------------------------------------------------------------------------
## Observations                      129        129         129         129    
## R2                               0.871      0.798       0.893       0.839   
## Adjusted R2                      0.870      0.796       0.892       0.837   
## Residual Std. Error (df = 127)   29.106     36.459      26.536      32.564  
## F Statistic (df = 1; 127)      858.827*** 501.275*** 1,058.984*** 660.555***
## ============================================================================
## Note:                                            *p<0.1; **p<0.05; ***p<0.01

Variability in \(R^2\): .798 for sensor 40 vs .892 for sensor 5. Also Sensor 40 “under-estimates” PAR the most and sensor 5 “over-estimates” PAR the most.

Effect of Depth 

p <- ggplot(hobo_licor,aes(x=Lux,y=PAR)) +
  facet_wrap(~Depth, scales = 'free') +
  geom_point(aes(col=Sample)) + 
  geom_smooth(formula = y ~ x, se=T, method = 'lm') + 
  theme_light() +
  theme(axis.text=element_text(size=8),
        axis.title=element_text(size=12))
ggplotly(p)

Figure 3 - PAR vs LUX and effect of depth on linear regression.

\(R^2\) of linear fit:

lm_fit_depth1 <- hobo_licor %>% 
  dplyr::filter(Depth == '1') %>% 
  lm(PAR ~ Lux, data = .) 

lm_fit_depth2 <- hobo_licor %>% 
  dplyr::filter(Depth == '2') %>% 
  lm(PAR ~ Lux, data = .) 

lm_fit_depth3 <- hobo_licor %>% 
  dplyr::filter(Depth == '3') %>% 
  lm(PAR ~ Lux, data = .) 

lm_fit_depth4 <- hobo_licor %>% 
  dplyr::filter(Depth == '4') %>% 
  lm(PAR ~ Lux, data = .) 

lm_fit_depth5 <- hobo_licor %>% 
  dplyr::filter(Depth == '5') %>% 
  lm(PAR ~ Lux, data = .) 

lm_fit_depthascent <- hobo_licor %>% 
  dplyr::filter(Depth == 'ascent') %>% 
  lm(PAR ~ Lux, data = .) 

stargazer(lm_fit_depth1,lm_fit_depth2,lm_fit_depth3,lm_fit_depth4,lm_fit_depth5,lm_fit_depthascent, title="Regression Results for Depth 1m (1), 2m (2), 3m (3), 4m (4), 5m (5) and ascent (6).",type = "text", style = 'default')#"html"
## 
## Regression Results for Depth 1m (1), 2m (2), 3m (3), 4m (4), 5m (5) and ascent (6).
## ===================================================================================================================================================================
##                                                                                   Dependent variable:                                                              
##                     -----------------------------------------------------------------------------------------------------------------------------------------------
##                                                                                           PAR                                                                      
##                               (1)                     (2)                     (3)                     (4)                     (5)                     (6)          
## -------------------------------------------------------------------------------------------------------------------------------------------------------------------
## Lux                        0.023***                0.022***                0.030***                0.024***                 0.026***                0.034***       
##                             (0.002)                 (0.002)                 (0.002)                 (0.002)                 (0.001)                 (0.007)        
##                                                                                                                                                                    
## Constant                  123.420***              114.616***               67.601***               81.741***               60.567***               85.694***       
##                            (11.897)                (13.775)                (11.400)                 (6.958)                 (3.003)                 (26.357)       
##                                                                                                                                                                    
## -------------------------------------------------------------------------------------------------------------------------------------------------------------------
## Observations                  84                      60                      48                      48                      244                      32          
## R2                           0.736                   0.730                   0.794                   0.808                   0.771                   0.427         
## Adjusted R2                  0.733                   0.725                   0.789                   0.804                   0.770                   0.408         
## Residual Std. Error    42.998 (df = 82)        45.064 (df = 58)        23.672 (df = 46)        14.624 (df = 46)        14.290 (df = 242)        51.643 (df = 30)   
## F Statistic         228.329*** (df = 1; 82) 156.511*** (df = 1; 58) 177.102*** (df = 1; 46) 194.013*** (df = 1; 46) 813.218*** (df = 1; 242) 22.326*** (df = 1; 30)
## ===================================================================================================================================================================
## Note:                                                                                                                                   *p<0.1; **p<0.05; ***p<0.01

Highest \(R^2\) at depth 4m, lowest for ascent (makes sense). Presence of “clusters” indicative of 2 different time of the day (morning vs afternoon).

Effect of Depth and Hobo sensor 

p <- ggplot(hobo_licor,aes(x=Lux,y=PAR,color=Hobo,group=Hobo)) +
  facet_wrap(~Depth, scales = 'free') +
  geom_point() + 
  geom_smooth(formula = y ~ x, se=T, method = 'lm') + 
  theme_light() +
  theme(axis.text=element_text(size=8),
        axis.title=element_text(size=12),
        strip.text = element_text(size=12))
ggplotly(p)

Figure 4 - PAR vs LUX and effect of sensor type and depth on linear regression.

Limitation:

  • I have assume a linear relationships between LUX and PAR, but can Long et al. (2012) help find a more common/generic relationship between LUX and PAR?

Non-linear regression PAR vs LUX

In Long et al. (2012),they assume a saturating behavior of PAR values (LiCOR) for high LUX values (from HOBO). They use the equation: \(PAR_{LiCOR} = A_1 \exp(-HOBO/t_1) + y_0\).

model_long <- nls(PAR ~ A1*exp(-Lux/t1) + y0, data=hobo_licor, start = list(A1 = -20000, t1 = 20000, y0 = 2000))

summary(model_long)
## 
## Formula: PAR ~ A1 * exp(-Lux/t1) + y0
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## A1  -592.66      48.75 -12.156  < 2e-16 ***
## t1 12840.07    1907.28   6.732 4.49e-11 ***
## y0   618.75      54.58  11.337  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.75 on 513 degrees of freedom
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 2.799e-06

Coefficients are compared to Long et al. (2012):

Quite lower than Long et al. (2012), not same range of LUX might be why + data quantity.

p <- ggplot() +
  geom_point(data=hobo_licor,aes(x=Lux,y=PAR)) +
  geom_line(data=hobo_licor,aes(x=Lux,y=predict(model_long)),col='red',size=1) +
  geom_line(data=hobo_licor,aes(x=Lux,y=predict(lm_fit)),col='blue',size=1)
ggplotly(p)

Figure 5 - PAR vs LUX, linear regression (blue) and exponential regression following Long et al. (2012) (red).

Let’s compare the RMSE between models:

paste0('Long Expo Model RMSE =',rmse(hobo_licor$PAR,predict(model_long)))
## [1] "Long Expo Model RMSE =32.652806338231"
paste0('Linear Model RMSE =',rmse(hobo_licor$PAR,predict(lm_fit)))
## [1] "Linear Model RMSE =34.4007905909135"

And let’s compare the AIC between models:

AIC(model_long,lm_fit)
##            df      AIC
## model_long  4 5069.825
## lm_fit      3 5121.643

The exponential model seems better. But we porbably need more data at higher LUX.

Bibliography

Long, Matthew H, Jennie E Rheuban, Peter Berg, and Joseph C Zieman. 2012. “A Comparison and Correction of Light Intensity Loggers to Photosynthetically Active Radiation Sensors.” Limnology and Oceanography: Methods 10 (6): 416–24.