Jack’s Analysis Results

Script Objective: This script will analyze the metrics derived from our metrics_combined csv

Libraries

library(ggplot2)
library(ggpmisc)

Loading required package: ggpp

Registered S3 methods overwritten by 'ggpp':
  method                  from   
  heightDetails.titleGrob ggplot2
  widthDetails.titleGrob  ggplot2

Attaching package: 'ggpp'

The following object is masked from 'package:ggplot2':

    annotate

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(lme4)

Loading required package: Matrix

library(patchwork)

Globals

out = read.csv("S:/ursa/jones/NRD/metrics_combined_mc.csv")
vegTypes = c("BTMAPLE","DRYFIR","GAMBLOAK","GBRPJ","LODGEP","QASPEN","SPRFIR")
radii = c(1,2,3,4,5,6,7,8)
figurePath = "S:/ursa/jones/Thesis/Figures"
custom_labels = c(
  "GBRPJ" = "Piñon-Juniper",
  "LODGEP" = "Lodgepole Pine",
  "SPRFIR" = "Spruce Mixed Conifer",
  "DRYFIR" = "Dry Mixed Conifer",
  "GAMBLOAK" = "Gambel Oak",
  "QASPEN" = "Quaking Aspen",
  "BTMAPLE" = "Bigtooth Maple"
)

1. Compare MLS density to ALS and ULS NRD at each radius with all vegetation types

df.allVegR2s = data.frame(compType=character(), radius=double(), r2=double())
df.perVegR2s.mlsals = data.frame(compType=character(), vegType=character(), radius=double(), r2=double(), od=double(), ht=double())
df.perVegR2s.mlsuls = data.frame(compType=character(), vegType=character(), radius=double(), r2=double(), od=double(), ht=double())

for (radius in radii){
  temp = out[out$radius == radius, ]
  
  temp.mls = data.frame(temp$mls_vol_dens, group = "mls")
  temp.als = data.frame(temp$als_nrd_under, group = "als")
  temp.uls = data.frame(temp$uls_nrd_under, group = "uls")

  lm.temp.mlsals = lm(temp$mls_vol_dens ~ temp$als_nrd_under, data = temp)
  lm.temp.mlsuls = lm(temp$mls_vol_dens ~ temp$uls_nrd_under, data = temp)

  r2.temp.mlsals = summary(lm.temp.mlsals)$adj.r.squared
  r2.temp.mlsuls = summary(lm.temp.mlsuls)$adj.r.squared
  
  df.allVegR2s[nrow(df.allVegR2s)+1,] = c("MLS vs ALS", as.numeric(radius), as.numeric(r2.temp.mlsals))
  df.allVegR2s[nrow(df.allVegR2s)+1,] = c("MLS vs ULS", as.numeric(radius), as.numeric(r2.temp.mlsuls))

}

p = ggplot(out, aes(x=als_nrd_under, y=mls_vol_dens, color=als_nrd_under)) + 
    geom_point() + 
    geom_smooth(aes(group = 1), color="red", method = "lm", se = FALSE) +
  stat_poly_eq(
    aes(label = paste(..rr.label..),group=1),
    formula = y ~ x,
    parse = TRUE,
    size = 4
  ) +
    theme_minimal() +
    scale_x_continuous(limits = c(0, 1)) +
    scale_y_continuous(limits = c(0, .8)) +
    scale_color_viridis_c(option = "D") +
    facet_wrap(~ radius, 
           labeller = 
             as_labeller(function(x) paste("Radius:", x, "m")),
           scales="free") +
    labs(x = "ALS NRD", y = "MLS Density", color="ALS NRD Value")

plot(p)

Warning: The dot-dot notation (`..rr.label..`) was deprecated in ggplot2 3.4.0.
i Please use `after_stat(rr.label)` instead.

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 74 rows containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 74 rows containing non-finite outside the scale range
(`stat_poly_eq()`).

Warning: The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?

Warning: Removed 74 rows containing missing values or values outside the scale range
(`geom_point()`).

  p = ggplot(out, aes(x=uls_nrd_under, y=mls_vol_dens,color=uls_nrd_under)) + 
    geom_point() + 
    geom_smooth(aes(group = 1), color="red", method = "lm", se = FALSE) +
  stat_poly_eq(
    aes(label = paste(..rr.label..),group=1),
    formula = y ~ x,
    parse = TRUE,
    size = 4
  ) +
    facet_wrap(~ radius, 
           labeller = 
             as_labeller(function(x) paste("Radius:", x, "m")),
           scales="free") +
    theme_minimal() +
    scale_x_continuous(limits = c(0, 1)) +
    scale_y_continuous(limits = c(0, .8)) +
    scale_color_viridis_c(option = "D") +
    labs(x = "ULS NRD", y = "MLS Density", color="ULS NRD Value")

  plot(p)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 41 rows containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 41 rows containing non-finite outside the scale range
(`stat_poly_eq()`).

Warning: The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?

Warning: Removed 41 rows containing missing values or values outside the scale range
(`geom_point()`).

2. Compare R^2’s between ALS and ULS at each radius

p = ggplot(df.allVegR2s, aes(x = as.numeric(radius), y = as.numeric(r2), color = factor(compType))) +
  geom_point(size = 2, shape = 21, stroke = 1) +
  labs(x = "Radius", y = "R2", color = "Comparison Type") +
  theme_bw() +
  theme(panel.border = element_rect(color = "black", fill = NA, size = 1))

Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
i Please use the `linewidth` argument instead.

p

3. Radius 4 seems to have the best combination of ULS and ALS results

Radius 4 seems to have the best combination of ULS and ALS results based on 2.

4. Assess effects of OD and OH on model performance

for (vegType in vegTypes){
  temp = out[out$veg_type == vegType & out$radius == 4, ]

  lm.temp.mlsals = lm(mls_vol_dens ~ als_nrd_under + als_ord_over + als_ht_over, data = temp)
  lm.temp.mlsuls = lm(mls_vol_dens ~ uls_nrd_under + uls_ord_over + uls_ht_over, data = temp)

  df.perVegR2s.mlsals[nrow(df.perVegR2s.mlsals)+1,] = c("MLS vs ALS", vegType, radius, summary(lm.temp.mlsals)$adj.r.squared, mean(temp$als_ord_over), quantile(temp$als_ht_over, 0.9))
  df.perVegR2s.mlsuls[nrow(df.perVegR2s.mlsuls)+1,] = c("MLS vs ULS", vegType, radius, summary(lm.temp.mlsuls)$adj.r.squared, mean(temp$uls_ord_over), quantile(temp$uls_ht_over, 0.9))
}

p1 = ggplot(df.perVegR2s.mlsals, aes(x = as.numeric(od), y = as.numeric(r2), color = factor(vegType))) +
geom_point(size = 2, shape = 21, stroke = 1) +
    geom_smooth(aes(group = 1), color="red", method = "lm", se = FALSE) +
  stat_poly_eq(
    aes(label = paste(..rr.label..),group=1),
    formula = y ~ x,
    parse = TRUE,
    size = 4,
    label.x=0.05,
    label.y=0.05
  ) +
labs(x = "ALS OD", y = "R²", color = "Veg Type") +
theme_bw() +
theme(panel.border = element_rect(color = "black", fill = NA, size = 1))

p2 = ggplot(df.perVegR2s.mlsals, aes(x = as.numeric(ht), y = as.numeric(r2), color = factor(vegType))) +
geom_point(size = 2, shape = 21, stroke = 1) +
    geom_smooth(aes(group = 1), color="red", method = "lm", se = FALSE) +
  stat_poly_eq(
    aes(label = paste(..rr.label..),group=1),
    formula = y ~ x,
    parse = TRUE,
    size = 4,
    label.x=0.05,
    label.y=0.05
  ) +
labs(x = "ALS HT", y = "R²", color = "Veg Type") +
theme_bw() +
theme(panel.border = element_rect(color = "black", fill = NA, size = 1))


p3 = ggplot(df.perVegR2s.mlsuls, aes(x = as.numeric(od), y = as.numeric(r2), color = factor(vegType))) +
geom_point(size = 2, shape = 21, stroke = 1) +
    geom_smooth(aes(group = 1), color="red", method = "lm", se = FALSE) +
  stat_poly_eq(
    aes(label = paste(..rr.label..),group=1),
    formula = y ~ x,
    parse = TRUE,
    size = 4,
    label.x=0.05,
    label.y=0.05
  ) +
  labs(x = "ULS OD", y = "R²", color = "Veg Type") +
  theme_bw() +
  theme(panel.border = element_rect(color = "black", fill = NA, size = 1))


p4 = ggplot(df.perVegR2s.mlsuls, aes(x = as.numeric(ht), y = as.numeric(r2), color = factor(vegType))) +
  geom_point(size = 2, shape = 21, stroke = 1) +
  geom_smooth(aes(group = 1), color="red", method = "lm", se = FALSE) +
  stat_poly_eq(
    aes(label = paste(..rr.label..),group=1),
    formula = y ~ x,
    parse = TRUE,
    size = 4,
    label.x=0.05,
    label.y=0.05
  ) +
  labs(x = "ULS HT", y = "R²", color = "Veg Type") +
  theme_bw() +
  theme(panel.border = element_rect(color = "black", fill = NA, size = 1))

(p1 | p2) /
(p3 | p4)

`geom_smooth()` using formula = 'y ~ x'

Warning: The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?

`geom_smooth()` using formula = 'y ~ x'

Warning: The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?

`geom_smooth()` using formula = 'y ~ x'

Warning: The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?

`geom_smooth()` using formula = 'y ~ x'

Warning: The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?

5. Build LMs to predict MLS density with ALS and ULS NRD per vegetation type

p = ggplot(out, aes(x=as.numeric(als_nrd_under), y=as.numeric(mls_vol_dens),color=als_nrd_under)) + 
    geom_point() + 
    geom_smooth(aes(group = 1), color="red", method = "lm", se = FALSE) +
  stat_poly_eq(
    aes(label = paste(..rr.label..),group=1),
    formula = y ~ x,
    parse = TRUE,
    size = 4
  ) +
    theme_minimal() +
    scale_x_continuous(limits = c(0, 1)) +
    scale_y_continuous(limits = c(0, .8)) +
    scale_color_viridis_c(option = "D") +
    facet_wrap(~ veg_type,
           scales="free",
           labeller = labeller(veg_type = custom_labels)) +
    labs(x = "ALS NRD", y = "MLS Density", color="ALS NRD Value") +
    theme(panel.border = element_rect(color = "black", fill = NA, size = 1))

plot(p)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 74 rows containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 74 rows containing non-finite outside the scale range
(`stat_poly_eq()`).

Warning: The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?

Warning: Removed 74 rows containing missing values or values outside the scale range
(`geom_point()`).

p = ggplot(out, aes(x=as.numeric(uls_nrd_under), y=as.numeric(mls_vol_dens),color=uls_nrd_under)) + 
    geom_point() + 
    geom_smooth(aes(group = 1), color="red", method = "lm", se = FALSE) +
  stat_poly_eq(
    aes(label = paste(..rr.label..),group=1),
    formula = y ~ x,
    parse = TRUE,
    size = 4
  ) +
    theme_minimal() +
    scale_x_continuous(limits = c(0, 1)) +
    scale_y_continuous(limits = c(0, .8)) +
    scale_color_viridis_c(option = "D") +
    facet_wrap(~ veg_type,
           scales="free",
           labeller = labeller(veg_type = custom_labels)) +
    labs(x = "ULS NRD", y = "MLS Density", color="ULS NRD Value") +
    theme(panel.border = element_rect(color = "black", fill = NA, size = 1))

plot(p)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 41 rows containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 41 rows containing non-finite outside the scale range
(`stat_poly_eq()`).

Warning: The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?
The following aesthetics were dropped during statistical transformation:
colour.
i This can happen when ggplot fails to infer the correct grouping structure in
  the data.
i Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?

Warning: Removed 41 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_smooth()`).

6.1 Build LMMs with vegType as a random effect

{# {r} # LMM.mlsals.1 = lmer(mls_vol_dens ~ als_nrd_under + (1 | veg_type), data = out) # summary(LMM.mlsals.1) # # LMM.mlsuls.1 = lmer(mls_vol_dens ~ uls_nrd_under + (1 | veg_type), data = out) # summary(LMM.mlsuls.1) # # LMM.mlsals = lmer(mls_vol_dens ~ (als_nrd_under | veg_type), data = out) # summary(LMM.mlsals) # # LMM.mlsuls = lmer(mls_vol_dens ~ (uls_nrd_under | veg_type), data = out) # summary(LMM.mlsuls) # # LMM.mlsals = lmer(mls_vol_dens ~ als_nrd_under + (als_nrd_under | veg_type), data = out) # summary(LMM.mlsals) # # LMM.mlsuls = lmer(mls_vol_dens ~ uls_nrd_under + (uls_nrd_under | veg_type), data = out) # summary(LMM.mlsuls)

6.2 lmer(mls_vol_dens ~ (als_nrd_under | veg_type), data = df_train)

# read in the data
df = out

# get list of plots
plot_ids = unique(df$plot)

# loop through them
for (p in plot_ids){

  # split data into training and validation
  df_train = df[df$plot_id != p & df$radius==4, ] # all plots *except* the one left out
  df_valid = df[df$plot_id == p & df$radius==4,] # just data from the one plot left out

  # build a model with the training data
  mod.mlsals = lmer(mls_vol_dens ~ (als_nrd_under | veg_type), data = df_train)
  
  mod.mlsuls = lmer(mls_vol_dens ~ (uls_nrd_under | veg_type), data = df_train)

  # make predictions on the validation data
  pred.mlsals = predict(mod.mlsals, df_valid)
  pred.mlsuls = predict(mod.mlsuls, df_valid)
  
  # create data.frame of predictions and observations for this "fold"
  df_pred_obs_fold.mlsals = data.frame(pred = pred.mlsals, 
                                 obs = df_valid$mls_vol_dens)

  df_pred_obs_fold.mlsuls = data.frame(pred = pred.mlsuls, 
                                 obs = df_valid$mls_vol_dens)
  # compile data from all folds
  if (p == plot_ids[1]){

    # on the first iteration, just copy the data.frame
    df_pred_obs.mlsals = df_pred_obs_fold.mlsals
    df_pred_obs.mlsuls = df_pred_obs_fold.mlsuls
    
  } else {

    # on subsequent iterations, combine the data by row
    df_pred_obs.mlsals = rbind(df_pred_obs.mlsals, df_pred_obs_fold.mlsals)
    df_pred_obs.mlsuls = rbind(df_pred_obs.mlsuls, df_pred_obs_fold.mlsuls)

  }

}

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00207738 (tol = 0.002, component 1)

# get r-squared and rmse
r2.mlsals = 1 - sum((df_pred_obs.mlsals$obs - df_pred_obs.mlsals$pred)^2) / sum((df_pred_obs.mlsals$obs - mean(df_pred_obs.mlsals$obs))^2)

rmse.mlsals = sqrt(mean((df_pred_obs.mlsals$pred - df_pred_obs.mlsals$obs)^2))

r2.mlsuls = 1 - sum((df_pred_obs.mlsuls$obs - df_pred_obs.mlsuls$pred)^2) / sum((df_pred_obs.mlsuls$obs - mean(df_pred_obs.mlsuls$obs))^2)

rmse.mlsuls = sqrt(mean((df_pred_obs.mlsuls$pred - df_pred_obs.mlsuls$obs)^2))

#plot out predicted vs observed
p = ggplot(data=df_pred_obs.mlsals,aes(x=obs,y=pred)) + 
    geom_point() + 
    geom_smooth(color="red", method = "lm", se = FALSE) +
  stat_poly_eq(
    aes(label = paste(..rr.label..),group=1),
    formula = y ~ x,
    parse = TRUE,
    size = 4
  ) +
    theme_minimal() +
    scale_x_continuous(limits = c(0, .45)) +
    scale_y_continuous(limits = c(0, .45)) +
    geom_abline(slope = 1, 
                intercept = 0, 
                linetype = "dashed", 
                color = "blue", 
                linewidth = 1) +
    annotate("text",
             x = 0.0, 
             y = 0.45,
             label = paste0("RMSE: ", round(rmse.mlsals, 3)),
             hjust = 0, vjust = 5, size = 4) +
    labs(x = "Observed", y = "Predicted") +
    theme(panel.border = element_rect(color = "black", fill = NA, size = 1))

plot(p)

`geom_smooth()` using formula = 'y ~ x'

p = ggplot(data=df_pred_obs.mlsuls,aes(x=obs,y=pred)) + 
    geom_point() + 
    geom_smooth( color="red", method = "lm", se = FALSE) +
  stat_poly_eq(
    aes(label = paste(..rr.label..),group=1),
    formula = y ~ x,
    parse = TRUE,
    size = 4
  ) +
    theme_minimal() +
    scale_x_continuous(limits = c(0, .45)) +
    scale_y_continuous(limits = c(0, .45)) +
    geom_abline(slope = 1, 
                intercept = 0, 
                linetype = "dashed", 
                color = "blue", 
                linewidth = 1) +
    annotate("text",
             x = 0.0, 
             y = 0.45,
             label = paste0("RMSE: ", round(rmse.mlsuls, 3)),
             hjust = 0, vjust = 5, size = 4) +
    labs(x = "Observed", y = "Predicted") +
    theme(panel.border = element_rect(color = "black", fill = NA, size = 1))

plot(p)

`geom_smooth()` using formula = 'y ~ x'

6.3 lmer(mls_vol_dens ~ als_nrd_under + (als_nrd_under | veg_type), data = df_train)

{# {r} # # read in the data # df = out # # # get list of plots # plot_ids = unique(df$plot) # # # loop through them # for (p in plot_ids){ # # # split data into training and validation # df_train = df[df$plot_id != p & df$radius==4, ] # all plots *except* the one left out # df_valid = df[df$plot_id == p & df$radius==4,] # just data from the one plot left out # # # build a model with the training data # mod.mlsals = lmer(mls_vol_dens ~ als_nrd_under + (als_nrd_under | veg_type), data = df_train) # # mod.mlsuls = lmer(mls_vol_dens ~ uls_nrd_under + (uls_nrd_under | veg_type), data = df_train) # # # make predictions on the validation data # pred.mlsals = predict(mod.mlsals, df_valid) # pred.mlsuls = predict(mod.mlsuls, df_valid) # # # create data.frame of predictions and observations for this "fold" # df_pred_obs_fold.mlsals = data.frame(pred = pred.mlsals, # obs = df_valid$mls_vol_dens) # # df_pred_obs_fold.mlsuls = data.frame(pred = pred.mlsuls, # obs = df_valid$mls_vol_dens) # # compile data from all folds # if (p == plot_ids[1]){ # # # on the first iteration, just copy the data.frame # df_pred_obs.mlsals = df_pred_obs_fold.mlsals # df_pred_obs.mlsuls = df_pred_obs_fold.mlsuls # # } else { # # # on subsequent iterations, combine the data by row # df_pred_obs.mlsals = rbind(df_pred_obs.mlsals, df_pred_obs_fold.mlsals) # df_pred_obs.mlsuls = rbind(df_pred_obs.mlsuls, df_pred_obs_fold.mlsuls) # # } # # } # # # get r-squared and rmse # r2.mlsals = 1 - sum((df_pred_obs.mlsals$obs - df_pred_obs.mlsals$pred)^2) / sum((df_pred_obs.mlsals$obs - mean(df_pred_obs.mlsals$obs))^2) # # rmse.mlsals = sqrt(mean((df_pred_obs.mlsals$pred - df_pred_obs.mlsals$obs)^2)) # # r2.mlsuls = 1 - sum((df_pred_obs.mlsuls$obs - df_pred_obs.mlsuls$pred)^2) / sum((df_pred_obs.mlsuls$obs - mean(df_pred_obs.mlsuls$obs))^2) # # rmse.mlsuls = sqrt(mean((df_pred_obs.mlsuls$pred - df_pred_obs.mlsuls$obs)^2)) # # #plot out pred vs observed # p = ggplot(data=df_pred_obs.mlsals,aes(x=pred,y=obs)) + # geom_point() + # geom_smooth(color="red", method = "lm", se = FALSE) + # stat_poly_eq( # aes(label = paste(..rr.label..),group=1), # formula = y ~ x, # parse = TRUE, # size = 4 # ) + # theme_minimal() + # annotate("text", # x = min(df_pred_obs.mlsals$pred, na.rm = TRUE), # y = max(df_pred_obs.mlsals$obs, na.rm = TRUE), # label = paste0("RMSE: ", round(rmse.mlsals, 3)), # hjust = 0, vjust = 5, size = 4) + # labs(title = "MLS vs ALS Pred vs Obs LMM", x = "Predicted", y = "Observed") + # theme(panel.border = element_rect(color = "black", fill = NA, size = 1)) # # plot(p) # # p = ggplot(data=df_pred_obs.mlsuls,aes(x=pred,y=obs)) + # geom_point() + # geom_smooth( color="red", method = "lm", se = FALSE) + # stat_poly_eq( # aes(label = paste(..rr.label..),group=1), # formula = y ~ x, # parse = TRUE, # size = 4 # ) + # theme_minimal() + # annotate("text", # x = min(df_pred_obs.mlsuls$pred, na.rm = TRUE), # y = max(df_pred_obs.mlsuls$obs, na.rm = TRUE), # label = paste0("RMSE: ", round(rmse.mlsuls, 3)), # hjust = 0, vjust = 5, size = 4) + # labs(title = "MLS vs ULS Pred vs Obs LMM", x = "Predicted", y = "Observed") + # theme(panel.border = element_rect(color = "black", fill = NA, size = 1)) # # plot(p)