Yield curves from Alberta PSP data
Set up
Load up necessary libraries
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.2 --
## v ggplot2 3.4.0 v purrr 1.0.0
## v tibble 3.1.8 v dplyr 1.0.10
## v tidyr 1.2.1 v stringr 1.5.0
## v readr 2.1.3 v forcats 0.5.2
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()Read in PSP data, and examine its structure. There is an issue with “opening number” (columns 13 and 63), but we will not use that field in any case.
#setwd("G:/My Drive/Teaching/RENR431/AlbertaPSP_working/src")
setwd("C:/Users/gwa/Desktop/AlbertaPSP_working/src")
#setwd("/home/gwa/PSP/data")
pspCompiled <- read_csv(file = '../data/compiled.csv')## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
## dat <- vroom(...)
## problems(dat)## Rows: 193500 Columns: 119
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (21): company, fmu, fma, aspect, utm_zone, datum, natural_subregion, eco...
## dbl (91): company_plot_number, company_stand_number, establishment_year, est...
## lgl (7): sampling_unit_number, plot_comment, su, regen_decht, th_n_m, th_m,...
##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.str(pspCompiled)## spc_tbl_ [193,500 x 119] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ company : chr [1:193500] "GOA" "GOA" "GOA" "GOA" ...
## $ company_plot_number : num [1:193500] 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 ...
## $ company_stand_number: num [1:193500] 1 1 1 1 1 1 1 1 1 1 ...
## $ establishment_year : num [1:193500] 1960 1960 1960 1960 1960 1960 1960 1960 1960 1960 ...
## $ establishment_month : num [1:193500] 7 7 7 7 7 7 7 7 7 7 ...
## $ establishment_day : num [1:193500] 23 23 23 23 23 23 23 23 23 23 ...
## $ fmu : chr [1:193500] "W14" "W14" "W14" "W14" ...
## $ fma : chr [1:193500] "FMA7500020" "FMA7500020" "FMA7500020" "FMA7500020" ...
## $ ats_township : num [1:193500] 62 62 62 62 62 62 62 62 62 62 ...
## $ ats_range : num [1:193500] 19 19 19 19 19 19 19 19 19 19 ...
## $ ats_meridian : num [1:193500] 5 5 5 5 5 5 5 5 5 5 ...
## $ ats_section : num [1:193500] 31 31 31 31 31 31 31 31 31 31 ...
## $ opening_number : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ sampling_unit_number: logi [1:193500] NA NA NA NA NA NA ...
## $ topographic_position: num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ elevation : num [1:193500] 811 811 811 811 811 811 811 811 811 811 ...
## $ slope : num [1:193500] 2 2 2 2 2 2 2 2 2 2 ...
## $ aspect : chr [1:193500] "SW" "SW" "SW" "SW" ...
## $ x_coord : num [1:193500] 510771 510771 510771 510771 510771 ...
## $ y_coord : num [1:193500] 6028487 6028487 6028487 6028487 6028487 ...
## $ utm_zone : chr [1:193500] "UTM11" "UTM11" "UTM11" "UTM11" ...
## $ datum : chr [1:193500] "NAD83" "NAD83" "NAD83" "NAD83" ...
## $ latitude : num [1:193500] 54.4 54.4 54.4 54.4 54.4 ...
## $ longitude : num [1:193500] -117 -117 -117 -117 -117 ...
## $ natural_subregion : chr [1:193500] "CM" "CM" "CM" "CM" ...
## $ ecosite_guide : chr [1:193500] "N" "N" "N" "N" ...
## $ ecosite : chr [1:193500] NA NA NA NA ...
## $ ecosite_phase : chr [1:193500] NA NA NA NA ...
## $ plot_comment : logi [1:193500] NA NA NA NA NA NA ...
## $ shared : num [1:193500] 1 1 1 1 1 1 1 1 1 1 ...
## $ strata : chr [1:193500] "HwSw" "HwSw" "HwSw" "HwSw" ...
## $ substrata : chr [1:193500] "Pb" "Pb" "Pb" "Pb" ...
## $ est_yr : num [1:193500] 1960 1960 1960 1960 1960 1960 1960 1960 1960 1960 ...
## $ mmt_num : num [1:193500] 1 1 1 1 1 1 1 1 1 1 ...
## $ mmt_yr : num [1:193500] 1960 1960 1960 1960 1960 1960 1960 1960 1960 1960 ...
## $ mmt_mo : num [1:193500] 7 7 7 7 7 7 7 7 7 7 ...
## $ age : num [1:193500] 81 81 81 81 81 81 81 81 81 81 ...
## $ stand_type : num [1:193500] 1 1 1 1 1 1 1 1 1 1 ...
## $ scale : chr [1:193500] "con_dec" "con_dec" "con_dec" "species" ...
## $ species : chr [1:193500] "CO" "DE" "TO" "AW" ...
## $ sphRegenOnly : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ sphBH : num [1:193500] 691.7 306.3 998 108.7 39.5 ...
## $ sphD15 : num [1:193500] 691.7 306.3 998 108.7 39.5 ...
## $ sphD91 : num [1:193500] 652.2 306.3 958.5 108.7 39.5 ...
## $ ba : num [1:193500] 20.814 22.01 42.825 8.801 0.625 ...
## $ vol_0000 : num [1:193500] 167.04 190.07 357.12 86.6 3.55 ...
## $ vol_1307 : num [1:193500] 153.01 181.23 334.24 83.33 2.74 ...
## $ vol_1510 : num [1:193500] 141.99 177.84 319.83 82.32 2.02 ...
## $ biomass : num [1:193500] 79398 89934 169331 43342 2819 ...
## $ carbon : num [1:193500] 39699 44967 84666 21671 1409 ...
## $ topht_n : num [1:193500] NA NA NA 10 NA NA 8 NA NA 10 ...
## $ topht : num [1:193500] NA NA NA 23.9 13.7 ...
## $ topht_stat : chr [1:193500] NA NA NA "MP" ...
## $ age_n : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ age_tot : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ age_bh : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ age_max : num [1:193500] 81.4 81.4 81.4 81.4 81.4 ...
## $ age_harv : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ age_fire : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ age_avi : num [1:193500] 90 90 90 90 90 90 90 90 90 90 ...
## $ si_n : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ si_bh : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ opening : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ su : logi [1:193500] NA NA NA NA NA NA ...
## $ topo : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ elev : num [1:193500] 811 811 811 811 811 811 811 811 811 811 ...
## $ utm : chr [1:193500] "UTM11" "UTM11" "UTM11" "UTM11" ...
## $ natsub : chr [1:193500] "CM" "CM" "CM" "CM" ...
## $ ecoguide : chr [1:193500] "N" "N" "N" "N" ...
## $ eco : chr [1:193500] NA NA NA NA ...
## $ ecophs : chr [1:193500] NA NA NA NA ...
## $ stand_origin : chr [1:193500] "N" "N" "N" "N" ...
## $ plot_type : num [1:193500] 3 3 3 3 3 3 3 3 3 3 ...
## $ tree_area : num [1:193500] 1012 1012 1012 1012 1012 ...
## $ tree_tag : num [1:193500] 9.1 9.1 9.1 9.1 9.1 9.1 9.1 9.1 9.1 9.1 ...
## $ sap_area : num [1:193500] 1012 1012 1012 1012 1012 ...
## $ sap_dbh : num [1:193500] 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 ...
## $ sap_ht : num [1:193500] 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3 ...
## $ regen_area : num [1:193500] 0 0 0 0 0 0 0 0 0 0 ...
## $ regen_conht : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ regen_decht : logi [1:193500] NA NA NA NA NA NA ...
## $ regen_n : num [1:193500] 0 0 0 0 0 0 0 0 0 0 ...
## $ regendone : num [1:193500] 0 0 0 0 0 0 0 0 0 0 ...
## $ fire_orig : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ target : num [1:193500] NA NA NA 10 10 NA 10 NA NA 10 ...
## $ _FREQ_ : num [1:193500] NA NA NA 10 4 NA 6 NA NA 10 ...
## $ th_aa : num [1:193500] NA NA NA 23.9 13.7 ...
## $ th_n_a : num [1:193500] NA NA NA 10 4 NA 8 NA NA 10 ...
## $ th_mp : num [1:193500] NA NA NA 23.9 NA ...
## $ th_n_mp : num [1:193500] NA NA NA 10 NA NA 8 NA NA 10 ...
## $ th_mm : num [1:193500] NA NA NA 21.7 NA ...
## $ th_n_mm : num [1:193500] NA NA NA 3 NA NA NA NA NA 2 ...
## $ th_ms : num [1:193500] NA NA NA 21.7 NA ...
## $ th_n_ms : num [1:193500] NA NA NA 3 0 NA 0 NA NA 1 ...
## $ th_mpe : num [1:193500] NA NA NA 23.9 13.7 ...
## $ th_n_mpe : num [1:193500] NA NA NA 10 4 NA 6 NA NA 10 ...
## $ th_n_m : logi [1:193500] NA NA NA NA NA NA ...
## $ th_m : logi [1:193500] NA NA NA NA NA NA ...
## $ th_n_aa : logi [1:193500] NA NA NA NA NA NA ...
## [list output truncated]
## - attr(*, "spec")=
## .. cols(
## .. company = col_character(),
## .. company_plot_number = col_double(),
## .. company_stand_number = col_double(),
## .. establishment_year = col_double(),
## .. establishment_month = col_double(),
## .. establishment_day = col_double(),
## .. fmu = col_character(),
## .. fma = col_character(),
## .. ats_township = col_double(),
## .. ats_range = col_double(),
## .. ats_meridian = col_double(),
## .. ats_section = col_double(),
## .. opening_number = col_double(),
## .. sampling_unit_number = col_logical(),
## .. topographic_position = col_double(),
## .. elevation = col_double(),
## .. slope = col_double(),
## .. aspect = col_character(),
## .. x_coord = col_double(),
## .. y_coord = col_double(),
## .. utm_zone = col_character(),
## .. datum = col_character(),
## .. latitude = col_double(),
## .. longitude = col_double(),
## .. natural_subregion = col_character(),
## .. ecosite_guide = col_character(),
## .. ecosite = col_character(),
## .. ecosite_phase = col_character(),
## .. plot_comment = col_logical(),
## .. shared = col_double(),
## .. strata = col_character(),
## .. substrata = col_character(),
## .. est_yr = col_double(),
## .. mmt_num = col_double(),
## .. mmt_yr = col_double(),
## .. mmt_mo = col_double(),
## .. age = col_double(),
## .. stand_type = col_double(),
## .. scale = col_character(),
## .. species = col_character(),
## .. sphRegenOnly = col_double(),
## .. sphBH = col_double(),
## .. sphD15 = col_double(),
## .. sphD91 = col_double(),
## .. ba = col_double(),
## .. vol_0000 = col_double(),
## .. vol_1307 = col_double(),
## .. vol_1510 = col_double(),
## .. biomass = col_double(),
## .. carbon = col_double(),
## .. topht_n = col_double(),
## .. topht = col_double(),
## .. topht_stat = col_character(),
## .. age_n = col_double(),
## .. age_tot = col_double(),
## .. age_bh = col_double(),
## .. age_max = col_double(),
## .. age_harv = col_double(),
## .. age_fire = col_double(),
## .. age_avi = col_double(),
## .. si_n = col_double(),
## .. si_bh = col_double(),
## .. opening = col_double(),
## .. su = col_logical(),
## .. topo = col_double(),
## .. elev = col_double(),
## .. utm = col_character(),
## .. natsub = col_character(),
## .. ecoguide = col_character(),
## .. eco = col_character(),
## .. ecophs = col_character(),
## .. stand_origin = col_character(),
## .. plot_type = col_double(),
## .. tree_area = col_double(),
## .. tree_tag = col_double(),
## .. sap_area = col_double(),
## .. sap_dbh = col_double(),
## .. sap_ht = col_double(),
## .. regen_area = col_double(),
## .. regen_conht = col_double(),
## .. regen_decht = col_logical(),
## .. regen_n = col_double(),
## .. regendone = col_double(),
## .. fire_orig = col_double(),
## .. target = col_double(),
## .. `_FREQ_` = col_double(),
## .. th_aa = col_double(),
## .. th_n_a = col_double(),
## .. th_mp = col_double(),
## .. th_n_mp = col_double(),
## .. th_mm = col_double(),
## .. th_n_mm = col_double(),
## .. th_ms = col_double(),
## .. th_n_ms = col_double(),
## .. th_mpe = col_double(),
## .. th_n_mpe = col_double(),
## .. th_n_m = col_logical(),
## .. th_m = col_logical(),
## .. th_n_aa = col_logical(),
## .. totba = col_double(),
## .. totregden = col_double(),
## .. conba = col_double(),
## .. conregden = col_double(),
## .. decba = col_double(),
## .. decregden = col_double(),
## .. pinba = col_double(),
## .. pinregden = col_double(),
## .. sfba = col_double(),
## .. sfregden = col_double(),
## .. sbba = col_double(),
## .. sbregden = col_double(),
## .. pjba = col_double(),
## .. pjregden = col_double(),
## .. firba = col_double(),
## .. firregden = col_double(),
## .. bwba = col_double(),
## .. bwregden = col_double(),
## .. pbba = col_double(),
## .. pbregden = col_double()
## .. )
## - attr(*, "problems")=<externalptr>problems(pspCompiled)There’s a whole bunch of stuff we don’t need. Let’s keep only the relevant bits
psp <- pspCompiled |> select(company_plot_number,mmt_num,mmt_yr,species,
vol_1510,elevation,
latitude,longitude,
natsub,strata,scale,age,age_tot,
stand_origin)
str(psp)## tibble [193,500 x 14] (S3: tbl_df/tbl/data.frame)
## $ company_plot_number: num [1:193500] 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 ...
## $ mmt_num : num [1:193500] 1 1 1 1 1 1 1 1 1 1 ...
## $ mmt_yr : num [1:193500] 1960 1960 1960 1960 1960 1960 1960 1960 1960 1960 ...
## $ species : chr [1:193500] "CO" "DE" "TO" "AW" ...
## $ vol_1510 : num [1:193500] 141.99 177.84 319.83 82.32 2.02 ...
## $ elevation : num [1:193500] 811 811 811 811 811 811 811 811 811 811 ...
## $ latitude : num [1:193500] 54.4 54.4 54.4 54.4 54.4 ...
## $ longitude : num [1:193500] -117 -117 -117 -117 -117 ...
## $ natsub : chr [1:193500] "CM" "CM" "CM" "CM" ...
## $ strata : chr [1:193500] "HwSw" "HwSw" "HwSw" "HwSw" ...
## $ scale : chr [1:193500] "con_dec" "con_dec" "con_dec" "species" ...
## $ age : num [1:193500] 81 81 81 81 81 81 81 81 81 81 ...
## $ age_tot : num [1:193500] NA NA NA NA NA NA NA NA NA NA ...
## $ stand_origin : chr [1:193500] "N" "N" "N" "N" ...Filtering and mutating
Filter out the records we don’t need. We will be interested in the “con_dec” scale only for this exercise, and we’ll limit ourselves to the CM and DM natural subregions.
psp <- psp |>
filter(scale=="con_dec") |>
filter(natsub == "CM" | natsub == "DM")
dim(psp)## [1] 5517 14head(psp)tail(psp)The plots are broken down by subplot. We don’t want that.
psp <- psp |>
mutate(plotnum = floor(company_plot_number))|>
mutate(subplotnum = (company_plot_number - plotnum) * 10)
head(psp)tail(psp)let’s aggregate up to the plot level
psp_agg <- psp |>
arrange(plotnum,mmt_num,species) |>
group_by(plotnum, mmt_num, species) |>
summarise(latitude=mean(latitude),longitude=mean(longitude),
elevation=mean(elevation),vol_1510=mean(vol_1510),
age=mean(age),mmt_yr=mean(mmt_yr))## `summarise()` has grouped output by 'plotnum', 'mmt_num'. You can override
## using the `.groups` argument.head(psp_agg)dim(psp_agg)## [1] 3333 9pivot so that volumes by species group are all on one row.
psp_pivot <- psp_agg |>
pivot_wider(names_from = species, values_from = vol_1510)
head(psp_pivot)psp_pivotsummary(psp_pivot)## plotnum mmt_num latitude longitude
## Min. : 1 Min. :1.000 Min. :53.03 Min. :-119.9
## 1st Qu.: 376 1st Qu.:1.000 1st Qu.:55.07 1st Qu.:-117.4
## Median : 641 Median :3.000 Median :55.64 Median :-115.3
## Mean :3309 Mean :3.167 Mean :56.13 Mean :-115.0
## 3rd Qu.:8185 3rd Qu.:4.000 3rd Qu.:57.19 3rd Qu.:-113.1
## Max. :9033 Max. :9.000 Max. :59.19 Max. :-110.1
##
## elevation age mmt_yr CO
## Min. :302.2 Min. :-16.00 Min. :1960 Min. : 0.0000
## 1st Qu.:400.0 1st Qu.: 11.50 1st Qu.:1987 1st Qu.: 0.0000
## Median :608.0 Median : 34.00 Median :1993 Median : 0.9925
## Mean :567.7 Mean : 52.12 Mean :1993 Mean : 67.6428
## 3rd Qu.:644.0 3rd Qu.: 84.12 3rd Qu.:2001 3rd Qu.:117.4875
## Max. :898.0 Max. :211.50 Max. :2019 Max. :488.0192
## NA's :49 NA's :192
## DE TO
## Min. : 0.000 Min. : 0.00
## 1st Qu.: 0.000 1st Qu.: 0.00
## Median : 5.751 Median : 55.67
## Mean : 60.714 Mean :128.36
## 3rd Qu.: 92.172 3rd Qu.:247.85
## Max. :489.629 Max. :595.12
## Let’s recover the stratum information. It’s possible that subplots within a plot have different strata. The “strata” assigned to the plot is the one associated with the first subplot.
psp_strat <- psp |>
arrange(plotnum,mmt_num) |>
group_by(plotnum, mmt_num) |>
summarise(stratum=first(strata))## `summarise()` has grouped output by 'plotnum'. You can override using the
## `.groups` argument.head(psp_strat)dim(psp_strat)## [1] 1111 3add stratum back to the aggregated data
mergeCols <- c("plotnum","mmt_num")
psp_merged <- merge(psp_pivot,psp_strat, by = mergeCols)
head(psp_merged)get rid of records with negative ages and rename stuff
psp_merged <- psp_merged|>
filter(age >= 0) |>
rename(vol_tot=TO, vol_con=CO, vol_dec=DE)
summary(psp_merged)## plotnum mmt_num latitude longitude
## Min. : 1.0 Min. :1.000 Min. :53.03 Min. :-119.4
## 1st Qu.: 349.0 1st Qu.:2.000 1st Qu.:55.08 1st Qu.:-117.4
## Median : 681.5 Median :3.000 Median :55.64 Median :-115.3
## Mean :3629.7 Mean :3.412 Mean :56.21 Mean :-114.9
## 3rd Qu.:9006.0 3rd Qu.:5.000 3rd Qu.:57.82 3rd Qu.:-113.0
## Max. :9033.0 Max. :9.000 Max. :59.19 Max. :-110.1
##
## elevation age mmt_yr vol_con
## Min. :303 Min. : 0.00 Min. :1960 Min. : 0.0
## 1st Qu.:381 1st Qu.: 12.00 1st Qu.:1987 1st Qu.: 0.0
## Median :608 Median : 35.00 Median :1992 Median : 0.0
## Mean :562 Mean : 53.26 Mean :1992 Mean : 64.7
## 3rd Qu.:640 3rd Qu.: 85.00 3rd Qu.:2001 3rd Qu.:102.8
## Max. :898 Max. :211.50 Max. :2019 Max. :488.0
## NA's :10
## vol_dec vol_tot stratum
## Min. : 0.000 Min. : 0.00 Length:902
## 1st Qu.: 0.000 1st Qu.: 0.00 Class :character
## Median : 1.803 Median : 16.83 Mode :character
## Mean : 50.393 Mean :115.09
## 3rd Qu.: 76.335 3rd Qu.:234.82
## Max. :489.629 Max. :595.12
## Analysis
Let’s find unique combos of plot and stratum
psp_merged |>
select(plotnum,stratum) |>
distinct(plotnum,stratum) |>
group_by(stratum) |>
summarise(count=n())We will use two yield strata, mixed (Sw, SwHw, HwSw, and Hw), and pine-black spruce (Pl, PlHw, and Sb).
Graphical summary of permanent sample plot data
ggplot(psp_merged, aes(x = age, y = vol_tot, colour = stratum, group = plotnum)) +
geom_line() +
ylim(0,600) +
xlab("age (y)") +
ylab("volume (m³/ha)") +
ggtitle ("Total volume")
Coniferous volume
ggplot(psp_merged, aes(x = age, y = vol_con, colour = stratum, group = plotnum)) +
geom_line() +
ylim(0,600) +
xlab("age (y)") +
ylab("volume (m³/ha)") +
ggtitle ("Coniferous volume")
Deciduous volume
ggplot(psp_merged, aes(x = age, y = vol_dec, colour = stratum, group = plotnum)) +
geom_line() +
ylim(0,600) +
xlab("age (y)") +
ylab("volume (m³/ha)") +
ggtitle ("Deciduous volume")
Based on the 3 plots above, it appears that a sigmoidal function of some sort would work well for total and coniferous yield. However, deciduous volumes appear to trend downward after 80-100 years.
Replot with aggregated strata.
psp_merged <- psp_merged |>
mutate(stratum2 = recode(stratum,
"Hw" = "Mixed",
"HwSw" = "Mixed",
"Sw" = "Mixed",
"SwHw" = "Mixed",
"Pl" = "PlSb",
"PlHw" = "PlSb",
"Sb" = "PlSb")) |>
arrange(stratum2,age)
ggplot(psp_merged, aes(x = age, y = vol_tot, colour = stratum2, group = plotnum)) +
geom_line() +
ylim(0,600)+
xlab("age (y)") +
ylab("volume (m³/ha)") +
ggtitle ("Total volume")
ggplot(psp_merged, aes(x = age, y = vol_con, colour = stratum2, group = plotnum)) +
geom_line() +
ylim(0,600)+
xlab("age (y)") +
ylab("volume (m³/ha)") +
ggtitle ("Coniferous volume")
ggplot(psp_merged, aes(x = age, y = vol_dec, colour = stratum2, group = plotnum)) +
geom_line() +
ylim(0,600)+
xlab("age (y)") +
ylab("volume (m³/ha)") +
ggtitle ("Deciduous volume")
Regression analysis
Use nonlinear regression to estimated coefficients for Chapman-Richards yield curve.
Mixed stratum
mixed <- psp_merged |>
filter(stratum2 == "Mixed")
y <- mixed$vol_tot
x <- mixed$age
m <- nls(y~b0*(1-exp(-b1*x))^b2,start=list(b0=300,b1=0.027,b2=4.0))
plot(x,y,type='p',ylim=c(0,500))
par(new=TRUE)
plot(x,predict(m),col="red",type='l',axes=FALSE,ylim=c(0,500))
summary(m)##
## Formula: y ~ b0 * (1 - exp(-b1 * x))^b2
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b0 3.464e+02 4.880e+00 70.982 < 2e-16 ***
## b1 5.558e-02 4.029e-03 13.796 < 2e-16 ***
## b2 1.730e+01 3.713e+00 4.661 3.7e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 54.65 on 789 degrees of freedom
##
## Number of iterations to convergence: 12
## Achieved convergence tolerance: 3.104e-06coef(m)## b0 b1 b2
## 346.39508943 0.05558123 17.30446221coeff_table <- tibble_row(stratum="mixed", type="tot",b0=coef(m)[1],b1=coef(m)[2],b2=coef(m)[3])mixed <- psp_merged |>
filter(stratum2 == "Mixed")
y <- mixed$vol_con
x <- mixed$age
m <- nls(y~b0*(1-exp(-b1*x))^b2,start=list(b0=300,b1=0.027,b2=4.0))
plot(x,y,type='p',ylim=c(0,500))
par(new=TRUE)
plot(x,predict(m),col="red",type='l',axes=FALSE,ylim=c(0,500))
summary(m)##
## Formula: y ~ b0 * (1 - exp(-b1 * x))^b2
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b0 2.801e+02 9.250e+00 30.284 <2e-16 ***
## b1 3.984e-02 4.145e-03 9.611 <2e-16 ***
## b2 2.535e+01 8.323e+00 3.046 0.0024 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 51.35 on 789 degrees of freedom
##
## Number of iterations to convergence: 10
## Achieved convergence tolerance: 2.441e-06coeff_table <- coeff_table |> add_row(stratum="mixed", type="con",b0=coef(m)[1],b1=coef(m)[2],b2=coef(m)[3])Pine/black spruce stratum
plsb <- psp_merged |>
filter(stratum2 == "PlSb")
y <- plsb$vol_tot
x <- plsb$age
m <- nls(y~300*(1-exp(-b1*x))^b2,start=list(b1=0.027,b2=4.0),control=nls.control(maxiter = 1000))
#m <- nls(y~b0*(1-exp(-b1*x))^b2,start=list(b0=450,b1=0.017702,b2=4.555790))
plot(x,y,type='p',ylim=c(0,500))
par(new=TRUE)
plot(x,predict(m),col="red",type='l',axes=FALSE,ylim=c(0,500))
summary(m)##
## Formula: y ~ 300 * (1 - exp(-b1 * x))^b2
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b1 0.02755 0.00446 6.178 1.17e-08 ***
## b2 7.06801 2.39056 2.957 0.00382 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 64.2 on 108 degrees of freedom
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 3.34e-06coeff_table <- coeff_table |> add_row(stratum="plsb", type="tot",b0=300,b1=coef(m)[1],b2=coef(m)[2])y <- plsb$vol_con
x <- plsb$age
m <- nls(y~200*(1-exp(-b1*x))^b2,start=list(b1=0.027,b2=4.0),control=nls.control(maxiter = 1000))
#m <- nls(y~200*(1-exp(-b1*x))^b2,start=list(b1=0.027,b2=4.0),control=nls.control(maxiter = 1000))
plot(x,y,type='p',ylim=c(0,500))
par(new=TRUE)
plot(x,predict(m),col="red",type='l',axes=FALSE,ylim=c(0,500))
summary(m)##
## Formula: y ~ 200 * (1 - exp(-b1 * x))^b2
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## b1 0.029837 0.005457 5.467 2.97e-07 ***
## b2 6.762636 2.615799 2.585 0.0111 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 50.82 on 108 degrees of freedom
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 8.844e-06coeff_table <- coeff_table |> add_row(stratum="plsb", type="con",b0=200,b1=coef(m)[1],b2=coef(m)[2])
coeff_table