COforest {cgam}R Documentation

Colorado Forest Data Set

Description

This data set contains 9167 records of different species of live trees for 345 sampled forested plots measured in 2015.

Usage

data("COforest")

Format

A data frame with 9167 observations on the following 19 variables.

PLT_CN

Unique identifier of plot

STATECD

State code using Bureau of Census Federal Information Processing Standards (FIPS)

COUNTYCD

County code (FIPS)

ELEV_PUBLIC

Elevation (ft) extracted spatially using LON_PUBLIC/LAT_PUBLIC

LON_PUBLIC

Fuzzed longitude in decimal degrees using NAD83 datum

LAT_PUBLIC

Fuzzed latitude in decimal degrees using NAD83 datum

ASPECT

a numeric vector

SLOPE

a numeric vector

SUBP

Subplot number

TREE

Tree number within subplot

STATUSCD

Tree status (0:no status; 1:live tree; 2:dead tree; 3:removed)

SPCD

Species code

DIA

Current diameter (in)

HT

Total height (ft): estimated when broken or missing

ACTUALHT

Actual height (ft): measured standing or down

HTCD

Height method code (1:measured; 2:actual measured-length estimated; 3:actual and length estimated; 4:modeled

TREECLCD

Tree class code (2:growing-stock; 3:rough cull; 4:rotten cull)

CR

Compacted crown ratio (percentage)

CCLCD

Crown class (1:open grown; 2:domimant; 3:co-dominant; 4:intermediate; 5:overtopped)

Source

It is provided by Forest Inventory Analysis (FIA) National Program.

References

X. Liao and M. Meyer (2019). Estimation and Inference in Mixed-Effect Regression Models using Shape Constraints, with Application to Tree Height Estimation. (to appear in Journal of the Royal Statistical Society. Series C: Applied Statistics)

Examples

## Not run: 
  library(dplyr)
  library(tidyr)
  data(COforest)

  #re-grouping classes of CCLCD:
  #combine dominant (2) and co-dominant (3)
  #combine intermediate (4) and overtopped (5)
  COforest = COforest 
    mutate(CCLCD = replace(CCLCD, CCLCD == 5, 4))

  #make a list of species, each element is a small data frame for one species
  species = COforest 

  #get the subset for quaking aspen, which is the 4th element in the species list
  sub = species$data[[4]]
  #for quaking aspen, there are only two crown classes: dominant/co-dominant
  #and intermediate/overtopped
  table(sub$CCLCD)
  #   3    4
  #1400  217
  #for quaking aspen, there are only two tree clases: growing-stock and rough cull
  table(sub$TREECLCD)
  #2    3
  #1591   26

  #fit the model
  ansc = cgamm(log(HT)~s.incr.conc(DIA)+factor(CCLCD)+factor(TREECLCD)
  +(1|PLT_CN), reml=TRUE, data=sub)

  #check which classes are significant
  summary(ansc)

  #fixed-effect 95
  newData = data.frame(DIA=sub$DIA,CCLCD=sub$CCLCD,TREECLCD=sub$TREECLCD)
  pfit = predict(ansc, newData,interval='confidence')
  lower = pfit$lower
  upper = pfit$upper
  #we need to use exp(muhat) later in the plot
  muhat = pfit$fit

  #get x and y
  x = sub$DIA
  y = sub$HT

  #get TREECLCD and CCLCD
  z1 = sub$TREECLCD
  z2 = sub$CCLCD

  #plot fixed-effect confidence intervals
  plot(x, y, xlab='DIA (m)', ylab='HT (m)', ylim=c(min(y),max(exp(upper))+10),type='n')
  lines(sort(x[z2==3&z1==2]), (exp(pfit$fit)[z2==3&z1==2])[order(x[z2==3&z1==2])],
  col='slategrey', lty=1, lwd=2)
  lines(sort(x[z2==3&z1==2]), (exp(pfit$lower)[z2==3&z1==2])[order(x[z2==3&z1==2])],
  col='slategrey', lty=1, lwd=2)
  lines(sort(x[z2==3&z1==2]), (exp(pfit$upper)[z2==3&z1==2])[order(x[z2==3&z1==2])],
  col='slategrey', lty=1, lwd=2)

  lines(sort(x[z2==4&z1==2]), (exp(pfit$fit)[z2==4&z1==2])[order(x[z2==4&z1==2])],
  col="blueviolet", lty=2, lwd=2)
  lines(sort(x[z2==4&z1==2]), (exp(pfit$lower)[z2==4&z1==2])[order(x[z2==4&z1==2])],
  col="blueviolet", lty=2, lwd=2)
  lines(sort(x[Cz2==4&z1==2]), (exp(pfit$upper)[Cz2==4&z1==2])[order(x[Cz2==4&z1==2])],
  col="blueviolet", lty=2, lwd=2)

  lines(sort(x[Cz2==3&z1==3]), (exp(pfit$fit)[Cz2==3&z1==3])[order(x[Cz2==3&z1==3])],
  col=3, lty=3, lwd=2)
  lines(sort(x[Cz2==3&z1==3]), (exp(pfit$lower)[Cz2==3&z1==3])[order(x[Cz2==3&z1==3])],
  col=3, lty=3, lwd=2)
  lines(sort(x[Cz2==3&z1==3]), (exp(pfit$upper)[Cz2==3&z1==3])[order(x[Cz2==3&z1==3])],
  col=3, lty=3, lwd=2)

  lines(sort(x[Cz2==4&z1==3]), (exp(pfit$fit)[Cz2==4&z1==3])[order(x[Cz2==4&z1==3])],
  col=2, lty=4, lwd=2)
  lines(sort(x[Cz2==4&z1==3]), (exp(pfit$lower)[Cz2==4&z1==3])[order(x[Cz2==4&z1==3])],
  col=2, lty=4, lwd=2)
  lines(sort(x[Cz2==4&z1==3]), (exp(pfit$upper)[Cz2==4&z1==3])[order(x[Cz2==4&z1==3])],
  col=2, lty=4, lwd=2)
  legend('bottomright', bty='n', cex=.9,c('dominant/co-dominant and growing stock',
  'intermediate/overtopped and growing stock','dominant/co-dominant and rough cull',
  'intermediate/overtopped and rough cull'), col=c('slategrey',"blueviolet",3,2),
  lty=c(1,2,3,4),lwd=c(2,2,2,2), pch=c(24,23,22,21))
  title('Quaking Aspen fits with 95

## End(Not run)
[Package cgam version 1.21 Index]