library(readxl)
library(lubridate)
library(tidyverse)
library(plotly)
library(stargazer)
library(Metrics)
## Warning: package 'Metrics' was built under R version 4.1.1

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 in Marlborough Sounds - data from Dan Crossett. 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

PAR vs LUX

All Depth and Sensors blend together 

hobo_licor <- read_excel('QCS_licor_hobo_calibration.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.023***          
##                               (0.001)          
##                                                
## Constant                    196.467***         
##                              (15.920)          
##                                                
## -----------------------------------------------
## Observations                    292            
## R2                             0.678           
## Adjusted R2                    0.677           
## Residual Std. Error     202.093 (df = 290)     
## F Statistic          610.844*** (df = 1; 290)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

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

\(PAR = 0.023*LUX + 196.467\)

Effect of Site 

p <- ggplot(hobo_licor,aes(x=Lux,y=PAR,col=Site,group=Site)) +
  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_siteDC <- hobo_licor %>% 
  dplyr::filter(Site == 'DC') %>% 
  lm(PAR ~ Lux, data = .) 

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

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

stargazer(lm_fit_siteDC,lm_fit_siteSB2,lm_fit_siteSC1, title="Regression Results for Sites DC (1), SB 2 (2) and SC 1 (3).",type = "text", style = 'default')#"html"
## 
## Regression Results for Sites DC (1), SB 2 (2) and SC 1 (3).
## ===========================================================================================
##                                               Dependent variable:                          
##                     -----------------------------------------------------------------------
##                                                       PAR                                  
##                               (1)                      (2)                    (3)          
## -------------------------------------------------------------------------------------------
## Lux                         0.021***                0.028***                0.019***       
##                             (0.002)                  (0.002)                (0.003)        
##                                                                                            
## Constant                   278.212***               58.328***              273.388***      
##                             (41.052)                (16.578)                (21.035)       
##                                                                                            
## -------------------------------------------------------------------------------------------
## Observations                  104                      94                      94          
## R2                           0.626                    0.789                  0.358         
## Adjusted R2                  0.622                    0.787                  0.351         
## Residual Std. Error    280.790 (df = 102)       111.229 (df = 92)      109.306 (df = 92)   
## F Statistic         170.419*** (df = 1; 102) 344.085*** (df = 1; 92) 51.338*** (df = 1; 92)
## ===========================================================================================
## Note:                                                           *p<0.1; **p<0.05; ***p<0.01

Variability in \(R^2\): 0.789 for site SB 2 vs 0.358 for site SC 1.

Effect of Hobo sensor 

p <- ggplot(hobo_licor,aes(x=Lux,y=PAR,col=as.character(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 3 - 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_sensor21 <- hobo_licor %>% 
  dplyr::filter(HOBO == 21) %>% 
  lm(PAR ~ Lux, data = .) 

stargazer(lm_fit_sensor21,lm_fit_sensor40, title="Regression Results for Sensors 21 (1), 40 (2).",type = "text", style = 'default')#"html"
## 
## Regression Results for Sensors 21 (1), 40 (2).
## ===========================================================
##                                    Dependent variable:     
##                                ----------------------------
##                                            PAR             
##                                     (1)            (2)     
## -----------------------------------------------------------
## Lux                               0.025***      0.022***   
##                                   (0.001)        (0.001)   
##                                                            
## Constant                         144.280***    236.249***  
##                                   (23.301)      (21.522)   
##                                                            
## -----------------------------------------------------------
## Observations                        146            146     
## R2                                 0.708          0.667    
## Adjusted R2                        0.706          0.665    
## Residual Std. Error (df = 144)    192.985        206.436   
## F Statistic (df = 1; 144)        348.612***    288.821***  
## ===========================================================
## Note:                           *p<0.1; **p<0.05; ***p<0.01

Variability in \(R^2\): .708 for sensor 21 vs .667 for sensor 40.

Effect of Depth 

p <- ggplot(hobo_licor,aes(x=Lux,y=PAR)) +
  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))
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.017***               0.014***               0.026***               0.025***                0.033***               0.002      
##                            (0.003)                (0.002)                (0.004)                (0.004)                 (0.002)               (0.021)     
##                                                                                                                                                           
## Constant                  486.388***             437.161***              113.277*              175.608***              75.124***             333.179**    
##                            (89.256)               (38.852)               (58.028)               (48.420)                (13.128)             (103.203)    
##                                                                                                                                                           
## ----------------------------------------------------------------------------------------------------------------------------------------------------------
## Observations                  40                     48                     26                     36                     136                    6        
## R2                          0.492                  0.644                  0.647                  0.593                   0.644                 0.002      
## Adjusted R2                 0.478                  0.636                  0.632                  0.581                   0.641                -0.247      
## Residual Std. Error   345.402 (df = 38)      167.527 (df = 46)      205.991 (df = 24)      203.369 (df = 34)       55.120 (df = 134)      54.301 (df = 4) 
## F Statistic         36.733*** (df = 1; 38) 83.089*** (df = 1; 46) 43.955*** (df = 1; 24) 49.544*** (df = 1; 34) 242.442*** (df = 1; 134) 0.010 (df = 1; 4)
## ==========================================================================================================================================================
## Note:                                                                                                                          *p<0.1; **p<0.05; ***p<0.01

Highest \(R^2\) at depth 5m, 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=as.character(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  -1725.8      145.6 -11.856  < 2e-16 ***
## t1  38386.4     5988.0   6.411 5.89e-10 ***
## y0   1791.5      157.0  11.411  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 180.8 on 289 degrees of freedom
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 1.907e-06

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

Quite lower than Long et al. (2012) except t1, 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 =179.868979906591"
paste0('Linear Model RMSE =',rmse(hobo_licor$PAR,predict(lm_fit)))
## [1] "Linear Model RMSE =201.400055220429"

And let’s compare the AIC between models:

AIC(model_long,lm_fit)
##            df      AIC
## model_long  4 3868.922
## lm_fit      3 3932.951

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.