This is the quantitative analysis of the 2013 National Woodland Owner (NWOS) Science Module on Landowner Values, fielded in Arkansas and Massachusetts. It is build around several questions intended to complement core questions on forest products and recreation, in order to be able to capture the Total Economic Value (TEV) of landowners’ wooded land. The paper will include some population-level estimates of key variables as well as these analyses.

First, we download the raw data from NWOS-DB:

library(nwos)
library(knitr)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## âś” dplyr     1.1.4     âś” readr     2.1.5
## âś” forcats   1.0.0     âś” stringr   1.5.1
## âś” ggplot2   3.5.1     âś” tibble    3.2.1
## âś” lubridate 1.9.4     âś” tidyr     1.3.1
## âś” purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## âś– dplyr::filter() masks stats::filter()
## âś– dplyr::lag()    masks stats::lag()
## â„ą Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(scales)
## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
library(DescTools)

citation(package="base") #citations
## To cite R in publications use:
## 
##   R Core Team (2024). _R: A Language and Environment for Statistical
##   Computing_. R Foundation for Statistical Computing, Vienna, Austria.
##   <https://www.R-project.org/>.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {R: A Language and Environment for Statistical Computing},
##     author = {{R Core Team}},
##     organization = {R Foundation for Statistical Computing},
##     address = {Vienna, Austria},
##     year = {2024},
##     url = {https://www.R-project.org/},
##   }
## 
## We have invested a lot of time and effort in creating R, please cite it
## when using it for data analysis. See also 'citation("pkgname")' for
## citing R packages.
R.Version()$version.string
## [1] "R version 4.4.2 (2024-10-31 ucrt)"
citation(package="DescTools")
## To cite package 'DescTools' in publications use:
## 
##   Signorell A (2025). _DescTools: Tools for Descriptive Statistics_. R
##   package version 0.99.59,
##   <https://CRAN.R-project.org/package=DescTools>.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {DescTools: Tools for Descriptive Statistics},
##     author = {Andri Signorell},
##     year = {2025},
##     note = {R package version 0.99.59},
##     url = {https://CRAN.R-project.org/package=DescTools},
##   }
#changing global settings
options(stringsAsFactors = FALSE)
options(scipen=999)

# quest <- get_nwos("2023","base intensified",c("25","5"),strict.intensification=T)
# plots <- get_nwos_plots_intensified("2023","base intensified",c("25","5"),strict.intensification=T)
# save.image("nwosdata.RData")
load("nwosdata.RData")

Sample size and cooperation rate (FFO only):

#response codes, family
rc <- plots@response_codes[plots@response_codes$OWNCD_NWOS=='45',]
table(rc$STATECD_NWOS,rc$RESPONSE_CAT) #response categories by state
##     
##        I  NC   P   R  UN
##   25 256  97  28 814 377
##   5  242  86  20 905 102
aggregate(RESPONSE_CAT~STATECD_NWOS,rc,"nwos_cooperation_rate") #cooperation rate by state
##   STATECD_NWOS RESPONSE_CAT
## 1           25    0.2331512
## 2            5    0.2073693

Next we download the standard estimates export. This includes all questions on the instrument, the core NWOS questions and the scimod-specific questions, as well as some custom question combinations:

ests <- readRDS("DATA/NWOS_2023_FFO_VALUES_ONE_PLUS.RDS")
ests$SE <- sqrt(ests$VARIANCE)

Analysis of the estimates

Next we look at some of the basic estimates for the SCIMOD questions.

#basic population summary

#total acres and ownerships
bps <- ests[ests$VARIABLE=='TOTAL'&ests$STATISTIC=='TOTAL',]
kable(bps[c('GEO_ABB','VARIABLE','UNITS','VALUE','SE')])
GEO_ABB VARIABLE UNITS VALUE SE
AR TOTAL ACRES 8516088.1 200280.52
AR TOTAL OWNERSHIPS 221766.6 39518.55
MA TOTAL ACRES 1336594.8 38161.89
MA TOTAL OWNERSHIPS 273600.1 35530.09 
#mean / median of AC_WOOD
ac <- ests[ests$VARIABLE=='AC_WOOD'&ests$STATISTIC %in% c('MEAN','MEDIAN'),]
kable(ac[c('GEO_ABB','VARIABLE','STATISTIC','VALUE','SE')])
GEO_ABB VARIABLE STATISTIC VALUE SE
AR AC_WOOD MEAN 38.401136 6.1051458
AR AC_WOOD MEDIAN 10.000000 3.0339993
MA AC_WOOD MEAN 4.885214 0.5574137
MA AC_WOOD MEDIAN 1.000000 0.3588703 
#mean / median of INC_WOOD
inc <- ests[ests$VARIABLE=='INC_WOOD'&ests$STATISTIC %in% c('MEAN','MEDIAN'),]
kable(inc[c('GEO_ABB','VARIABLE','STATISTIC','VALUE','SE')])
GEO_ABB VARIABLE STATISTIC VALUE SE
AR INC_WOOD MEAN 0.3579262 0.1093182
AR INC_WOOD MEDIAN 0.0000000 0.0000000
MA INC_WOOD MEAN 0.0887591 0.0515988
MA INC_WOOD MEDIAN 0.0000000 0.0000000 
#TEV questions

vsm <- ests[ests$STATISTIC=='PROPORTION'&
               !is.na(ests$UNITS)&
               ests$UNITS=='OWNERSHIPS',]

#identify variables to include
vsm.vars <- c('VSM_BEAUT','VSM_PEACE','VSM_LEARN','VSM_NATURE',
              'VSM_FUTURE','VSM_OTH_FUT','VSM_OTHERS','VSM_EXIST')

vsm <- vsm[vsm$VARIABLE %in% vsm.vars,]

#rename / label
vsm_arr <- c("Does not know",
                  "Does not enjoy",
                  "Slightly enjoys",
                  "Somewhat enjoys",
                  "Moderately enjoys",
                  "Greatly enjoys")
vsm$LABEL <- vsm_arr[match(vsm$LEVEL,c(9,1:5))]
vsm$LABEL <- factor(vsm$LABEL,levels=vsm_arr)

vsm$VARIABLE <- factor(vsm$VARIABLE,
                               levels = vsm.vars,
                               labels = c("Beauty", 
                                          "Peace",
                                          "Learning",
                                          "Nature",
                                          "Option",
                                          "Bequest",
                                          "Altruist",
                                          "Existence"))

vsm <- vsm[order(vsm$GEO_ABB,vsm$VARIABLE,vsm$LABEL),]
kable(vsm[c('GEO_ABB','VARIABLE','LABEL','VALUE','SE')])
GEO_ABB VARIABLE LABEL VALUE SE
AR Beauty Does not know 0.0038566 0.0025319
AR Beauty Does not enjoy 0.0000000 0.0000000
AR Beauty Slightly enjoys 0.0104330 0.0055973
AR Beauty Somewhat enjoys 0.1114421 0.0514300
AR Beauty Moderately enjoys 0.0730460 0.0325171
AR Beauty Greatly enjoys 0.8012224 0.0607140
AR Peace Does not know 0.0093221 0.0051025
AR Peace Does not enjoy 0.0038350 0.0022281
AR Peace Slightly enjoys 0.0013539 0.0019181
AR Peace Somewhat enjoys 0.0942308 0.0527285
AR Peace Moderately enjoys 0.0639208 0.0260151
AR Peace Greatly enjoys 0.8273375 0.0569546
AR Learning Does not know 0.1256846 0.0351536
AR Learning Does not enjoy 0.0039973 0.0023833
AR Learning Slightly enjoys 0.1312789 0.0561753
AR Learning Somewhat enjoys 0.2675275 0.0767592
AR Learning Moderately enjoys 0.1071971 0.0275796
AR Learning Greatly enjoys 0.3643146 0.0802040
AR Nature Does not know 0.0357189 0.0185461
AR Nature Does not enjoy 0.0000000 0.0000000
AR Nature Slightly enjoys 0.0205753 0.0173185
AR Nature Somewhat enjoys 0.0970760 0.0516733
AR Nature Moderately enjoys 0.1954079 0.0530361
AR Nature Greatly enjoys 0.6512219 0.0714780
AR Option Does not know 0.0363666 0.0181904
AR Option Does not enjoy 0.0002687 0.0004725
AR Option Slightly enjoys 0.0549640 0.0366706
AR Option Somewhat enjoys 0.1179974 0.0405036
AR Option Moderately enjoys 0.2774245 0.0668173
AR Option Greatly enjoys 0.5129788 0.0745758
AR Bequest Does not know 0.1253048 0.0373042
AR Bequest Does not enjoy 0.0322071 0.0108416
AR Bequest Slightly enjoys 0.1543655 0.0567448
AR Bequest Somewhat enjoys 0.0590725 0.0216981
AR Bequest Moderately enjoys 0.1935095 0.0480118
AR Bequest Greatly enjoys 0.4355405 0.0795333
AR Altruist Does not know 0.0795896 0.0223198
AR Altruist Does not enjoy 0.0243698 0.0094430
AR Altruist Slightly enjoys 0.1742469 0.0652620
AR Altruist Somewhat enjoys 0.1061283 0.0378021
AR Altruist Moderately enjoys 0.1435012 0.0396782
AR Altruist Greatly enjoys 0.4721641 0.0756313
AR Existence Does not know 0.0451882 0.0207983
AR Existence Does not enjoy 0.0134946 0.0099691
AR Existence Slightly enjoys 0.0120165 0.0050550
AR Existence Somewhat enjoys 0.1169648 0.0493586
AR Existence Moderately enjoys 0.1382622 0.0469237
AR Existence Greatly enjoys 0.6740735 0.0712552
MA Beauty Does not know 0.0004539 0.0005735
MA Beauty Does not enjoy 0.0000000 0.0000000
MA Beauty Slightly enjoys 0.0112130 0.0078259
MA Beauty Somewhat enjoys 0.0335728 0.0121316
MA Beauty Moderately enjoys 0.2138707 0.0470919
MA Beauty Greatly enjoys 0.7408895 0.0456434
MA Peace Does not know 0.0005242 0.0005871
MA Peace Does not enjoy 0.0004437 0.0007513
MA Peace Slightly enjoys 0.0047741 0.0032952
MA Peace Somewhat enjoys 0.0159703 0.0117418
MA Peace Moderately enjoys 0.1257025 0.0418215
MA Peace Greatly enjoys 0.8525852 0.0428096
MA Learning Does not know 0.1001741 0.0358570
MA Learning Does not enjoy 0.0265755 0.0136300
MA Learning Slightly enjoys 0.1383882 0.0433404
MA Learning Somewhat enjoys 0.3018827 0.0515490
MA Learning Moderately enjoys 0.2038093 0.0465155
MA Learning Greatly enjoys 0.2291701 0.0560809
MA Nature Does not know 0.0027673 0.0025314
MA Nature Does not enjoy 0.0002677 0.0002753
MA Nature Slightly enjoys 0.0037262 0.0025890
MA Nature Somewhat enjoys 0.0786395 0.0366988
MA Nature Moderately enjoys 0.1956800 0.0483007
MA Nature Greatly enjoys 0.7189194 0.0562962
MA Option Does not know 0.0703397 0.0259922
MA Option Does not enjoy 0.0128026 0.0133277
MA Option Slightly enjoys 0.0202820 0.0115440
MA Option Somewhat enjoys 0.1723581 0.0444463
MA Option Moderately enjoys 0.2548329 0.0515426
MA Option Greatly enjoys 0.4693848 0.0535427
MA Bequest Does not know 0.1381025 0.0393716
MA Bequest Does not enjoy 0.0706815 0.0287215
MA Bequest Slightly enjoys 0.1180198 0.0351428
MA Bequest Somewhat enjoys 0.2067635 0.0391816
MA Bequest Moderately enjoys 0.1790186 0.0441286
MA Bequest Greatly enjoys 0.2874141 0.0609888
MA Altruist Does not know 0.2000008 0.0475387
MA Altruist Does not enjoy 0.0143302 0.0096595
MA Altruist Slightly enjoys 0.0892977 0.0279809
MA Altruist Somewhat enjoys 0.2459706 0.0436508
MA Altruist Moderately enjoys 0.2008469 0.0445821
MA Altruist Greatly enjoys 0.2495538 0.0575531
MA Existence Does not know 0.0272577 0.0121108
MA Existence Does not enjoy 0.0000786 0.0001079
MA Existence Slightly enjoys 0.0039549 0.0025144
MA Existence Somewhat enjoys 0.0645054 0.0327544
MA Existence Moderately enjoys 0.2194982 0.0435209
MA Existence Greatly enjoys 0.6847052 0.0522118

We can also plot the full TEV, using a combination of the core NWOS variables pertaining to consumptive and non-consumptive direct values (binary variables) and those responding with the top two levels of the new scimod variables (i.e., Likert transformed as Box2 variables).

fig.tev <- ests[ests$STATISTIC=='PROPORTION'&
                !is.na(ests$UNITS)&
                ests$UNITS=='OWNERSHIPS',]
fig.tev$VALUE <- fig.tev$VALUE
fig.tev$SE <- sqrt(fig.tev$VARIANCE)

fig.tev <- fig.tev[fig.tev$GEO_ABB %in% c('AR','MA'),] #just two states
fig.tev$GEO_ABB <- factor(fig.tev$GEO_ABB)

# transition into plot

#identify variables to include and set labels
v.vars <- c('CUT','ACT_NTFP','ACT_GRAZE','FISHGAME',
            'RECNCONS','VSM_BEAUT_B2','VSM_PEACE_B2','VSM_LEARN_B2',
            'VSM_NATURE_B2','VSM_FUTURE_B2','VSM_OTH_FUT_B2','VSM_OTHERS_B2','VSM_EXIST_B2')
v.labels <- c('wood','NTFP','forage','fish/game',
              'recreation','beauty','peace and quiet','learning',
              'nature benefits','future benefits to self', 'future benefits to family','benefits to others','existence')

fig.tev <- fig.tev[fig.tev$VARIABLE %in% v.vars & fig.tev$LEVEL=='1',]
fig.tev$LABEL <- ordered(fig.tev$VARIABLE,labels=v.labels,levels=v.vars)

#categorize into TEV categories and set labels
vc <- c("DC","DNC","IN","OV","BV","AV","EV")
vc.labels <- c("direct (consumptive)","direct (non-consumptive)",
               "indirect","option","bequest","altruist","existence")

fig.tev$VC <- vc[1] #default is direct / consumptive
fig.tev$VC[fig.tev$VARIABLE %in% v.vars[5:8]] <- vc[2]
fig.tev$VC[fig.tev$VARIABLE %in% v.vars[9]] <- vc[3]
fig.tev$VC[fig.tev$VARIABLE %in% v.vars[10]] <- vc[4]
fig.tev$VC[fig.tev$VARIABLE %in% v.vars[11]] <- vc[5]
fig.tev$VC[fig.tev$VARIABLE %in% v.vars[12]] <- vc[6]
fig.tev$VC[fig.tev$VARIABLE %in% v.vars[13]] <- vc[7]

fig.tev$VC <- ordered(fig.tev$VC,labels=vc.labels,levels=vc)

#order variables

limits <- aes(ymax=(VALUE)+1*(SE),
              ymin=(VALUE)-1*(SE))

tev <- ggplot(data=fig.tev,aes(x=LABEL,y=VALUE,fill=GEO_ABB))+
  facet_grid(.~VC,scales="free",space="free",labeller = label_wrap_gen(width=10))+
  #coord_flip()+
  geom_bar(stat="identity",position="dodge")+
  geom_errorbar(limits,width=0,linewidth=0.5,position=position_dodge(width=0.90))+
  labs(x="value",y="percent of ownerships",fill="state")+
  scale_x_discrete(labels=wrap_format(7))+
  scale_y_continuous(label=percent,limits=c(0,1))+ #or whatever
  theme(text=element_text(size=12),
        panel.background=element_blank(),
        panel.grid=element_blank(),
        panel.border=element_rect(fill="transparent"),
        axis.ticks=element_line(color="black"),
        axis.text=element_text(color="black"))
print(tev)

svg("figure_3_TEV.svg",width=11,height=8)
print(tev)
dev.off()
## png 
##   2
fig.tev <- fig.tev[order(fig.tev$GEO_ABB,fig.tev$VC,fig.tev$LABEL),]
kable(fig.tev[c('GEO_ABB','VC','LABEL','VALUE','SE')])
GEO_ABB VC LABEL VALUE SE
AR direct (consumptive) wood 0.6113897 0.0777886
AR direct (consumptive) NTFP 0.1410705 0.0630447
AR direct (consumptive) forage 0.0722103 0.0255052
AR direct (consumptive) fish/game 0.4159042 0.0759199
AR direct (non-consumptive) recreation 0.4681346 0.0804709
AR direct (non-consumptive) beauty 0.8742683 0.0509118
AR direct (non-consumptive) peace and quiet 0.8912583 0.0520478
AR direct (non-consumptive) learning 0.4715117 0.0866777
AR indirect nature benefits 0.8466298 0.0545510
AR option future benefits to self 0.7904033 0.0529172
AR bequest future benefits to family 0.6290500 0.0656331
AR altruist benefits to others 0.6156653 0.0656625
AR existence existence 0.8123358 0.0529395
MA direct (consumptive) wood 0.6937492 0.0591277
MA direct (consumptive) NTFP 0.1176740 0.0381707
MA direct (consumptive) forage 0.0207695 0.0106205
MA direct (consumptive) fish/game 0.1109836 0.0381296
MA direct (non-consumptive) recreation 0.5042757 0.0621474
MA direct (non-consumptive) beauty 0.9547602 0.0156606
MA direct (non-consumptive) peace and quiet 0.9782877 0.0123223
MA direct (non-consumptive) learning 0.4329795 0.0587422
MA indirect nature benefits 0.9145993 0.0374560
MA option future benefits to self 0.7242177 0.0476260
MA bequest future benefits to family 0.4664327 0.0624920
MA altruist benefits to others 0.4504006 0.0654591
MA existence existence 0.9042034 0.0357694

Analysis of the raw data

First, we format the data:

wide <- nwos_wide(quest) #surveys, wide data.frame
wide[wide=='-1'] <- NA #-1 string equals null
wide <- wide[wide$OWNCD_NWOS=='45',] #just family

#does landowner hunt or fish? (same as construction as those in the estimates)
wide$FISHGAME <- ifelse(wide$REC_WHO_OWN=='1' &
                          (wide$REC_HOW_FISH=='1' |
                             wide$REC_HOW_HUNT=='1'),'1','0')

#does landowner engage in nonconsumptive recreation? (same construction as those in the estimates)
wide$RECNCONS <- ifelse(wide$REC_WHO_OWN=='1' &
                          (wide$REC_HOW_HIKE=='1' |
                             wide$REC_HOW_BIKE=='1' |
                             wide$REC_HOW_CAMP=='1' |
                             wide$REC_HOW_HORSE=='1' |
                             wide$REC_HOW_SKI=='1' |
                             wide$REC_HOW_OHV=='1'),'1','0')

Some chi-squared models to look at differences between states:

cs.vars <- c(c("CUT","ACT_NTFP","ACT_GRAZE",
               "FISHGAME","RECNCONS"),
             vsm.vars)
cs <- function(x){
  var <- wide[[x]] #choose variable
  xs <- chisq.test(wide$STATECD_NWOS,var,simulate.p.value=T)
  tt <- TschuprowT(wide$STATECD_NWOS,var,useNA="ifany")
  ret <- paste("X-squared = ",xs$statistic,", p-value = ",xs$p.value,", Tschuprow's t = ",tt,sep="")
  return(ret)
}
sapply(cs.vars,"cs")
##                                                                                               CUT 
##    "X-squared = 2.94149979045883, p-value = 0.105447276361819, Tschuprow's t = 0.064714288588512" 
##                                                                                          ACT_NTFP 
##    "X-squared = 2.4426768671195, p-value = 0.137431284357821, Tschuprow's t = 0.0609940512310795" 
##                                                                                         ACT_GRAZE 
## "X-squared = 17.8669388610005, p-value = 0.000499750124937531, Tschuprow's t = 0.160055541405527" 
##                                                                                          FISHGAME 
## "X-squared = 39.4868950742655, p-value = 0.000499750124937531, Tschuprow's t = 0.239135099284664" 
##                                                                                          RECNCONS 
##  "X-squared = 5.35691440482871, p-value = 0.0269865067466267, Tschuprow's t = 0.0935223087361076" 
##                                                                                         VSM_BEAUT 
##   "X-squared = 6.23866687569199, p-value = 0.193903048475762, Tschuprow's t = 0.0861050806931419" 
##                                                                                         VSM_PEACE 
##   "X-squared = 5.43476797536902, p-value = 0.385307346326837, Tschuprow's t = 0.0824627347502995" 
##                                                                                         VSM_LEARN 
##   "X-squared = 6.32017445781304, p-value = 0.272863568215892, Tschuprow's t = 0.0757589008090263" 
##                                                                                        VSM_NATURE 
##   "X-squared = 3.73939526704605, p-value = 0.600199900049975, Tschuprow's t = 0.0641297593105196" 
##                                                                                        VSM_FUTURE 
##  "X-squared = 17.8023538411781, p-value = 0.00349825087456272, Tschuprow's t = 0.124679970136576" 
##                                                                                       VSM_OTH_FUT 
##    "X-squared = 7.78962644288821, p-value = 0.15192403798101, Tschuprow's t = 0.0957740409883328" 
##                                                                                        VSM_OTHERS 
##   "X-squared = 3.88110113148967, p-value = 0.563718140929535, Tschuprow's t = 0.0779445320692226" 
##                                                                                         VSM_EXIST 
##   "X-squared = 8.25256995091671, p-value = 0.133433283358321, Tschuprow's t = 0.0901935541507158"

We also plot the raw sample of the new values questions…

#created stacked dataframe for plotting
long <- stack(wide[vsm.vars])
long$ind <- factor(long$ind,
                       levels = vsm.vars,
                       labels = c("Beauty", 
                                  "Peace",
                                  "Learning",
                                  "Nature",
                                  "Option",
                                  "Bequest",
                                  "Altruist",
                                  "Existence"))
long$values <- ordered(long$values,levels=c(9,1:5),labels=vsm_arr)
long$STATE <- rep(wide$STATECD_NWOS,length(vsm.vars)) #add state
long$STATE <- factor(long$STATE,levels=c(5,25),labels=c('AR','MA'))
long$val <- 1 #equal unit value
agg <- aggregate(val~values+ind+STATE,long,"length")

vsm.plot <- ggplot(agg,aes(x=ind,y=val,fill=values))+ #plots
  geom_bar(stat="identity")+
  coord_flip()+
  facet_grid(~STATE)+
  labs(x="values",y="frequency",fill="enjoyment")+
  theme(text=element_text(size=12),
        panel.background=element_blank(),
        panel.grid=element_blank(),
        panel.border=element_rect(fill="transparent"),
        axis.ticks=element_line(color="black"),
        axis.text=element_text(color="black"))
print(vsm.plot)

svg("figure_4_VSM.svg",width=11,height=8)
print(vsm.plot)
dev.off()
## png 
##   2