Marsh Billings Summary Statistics
Matthias Nevins
May 2, 2017
Read in data
# Read in CSV file
myData<-read.csv("MBR_Overstory_obs_wrking.csv",stringsAsFactors = F, header=T)
#str(myData)
#head(myData)
opar <- par(no.readonly=TRUE)
## Separate out observations on dead trees
mbdead <- myData[myData$tree_alive=="FALSE",]
## Now remove dead tree observations so the rest of your calculations are for live only.
myData2 <- myData[myData$tree_alive=="TRUE",]
#head(myData2)
attach(myData2)Transform data structure
# Convert to factors
myData2$STAND <- as.factor(myData2$STAND)
myData2$PLOT <- as.factor(myData2$PLOT)
myData2$tree_spp <- as.factor(myData2$tree_spp)
#str(myData2)Initial Observations
## Plot a histogram of tree DBH
hist(myData2$tree_dbh,xlab="Tree DBH - INCHES",nclass=40,col="blue") 
table(STAND) ## STAND
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## 41 24 27 56 9 33 44 34 16 9 19 5 38 16 17 11 63 66
## 19 20 21 22 23 24 25 26 27 28 30 31 32 33 34 35 37 38
## 16 30 64 19 11 35 27 29 15 16 40 60 19 45 53 19 41 46
## 39 40 41 42 43 44 45 46
## 180 36 15 34 15 47 53 72table(tree_spp)## tree_spp
## ACRU ACSA3 AMELA BEAL2 BELE BEPA FAGR FRAM2 JUCI LARIX OSVI PIAB
## 15 349 1 51 1 8 74 114 1 26 17 137
## PIRE PIRU PIST PISY POGR4 PRSE2 QURU ROPS TIAM TSCA ULAM
## 131 1 277 37 9 9 34 3 16 253 1table(PLOT)## PLOT
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## 14 8 15 14 6 7 9 10 6 5 13 13 12 12 13 8 10 9
## 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
## 6 3 7 13 6 11 13 3 7 7 15 5 3 7 5 6 6 11
## 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
## 13 7 10 10 15 11 9 14 8 16 14 14 6 2 8 12 23 15
## 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
## 10 7 14 7 13 6 6 6 6 4 6 13 8 9 10 13 23 7
## 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
## 16 9 7 19 11 6 7 8 13 18 3 13 14 13 15 10 10 7
## 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
## 17 10 16 12 9 14 8 16 1 3 1 13 4 10 14 6 7 8
## 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
## 16 10 8 7 12 17 16 20 13 20 11 13 17 4 5 5 2 9
## 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
## 8 13 13 8 4 4 8 1 18 10 12 12 13 9 2 8 23 16
## 145 146 147 148 150 151 152 153 154
## 12 13 2 7 15 17 12 16 23## Look at a table of the frequency of tree species
spptable <- table(tree_spp)
spptable## tree_spp
## ACRU ACSA3 AMELA BEAL2 BELE BEPA FAGR FRAM2 JUCI LARIX OSVI PIAB
## 15 349 1 51 1 8 74 114 1 26 17 137
## PIRE PIRU PIST PISY POGR4 PRSE2 QURU ROPS TIAM TSCA ULAM
## 131 1 277 37 9 9 34 3 16 253 1## assign that to a data.frame that is ordered by the least to most frequent species
spptable <- data.frame(spptable)
spptable <- spptable[order(spptable$Freq),]
spptable## tree_spp Freq
## 3 AMELA 1
## 5 BELE 1
## 9 JUCI 1
## 14 PIRU 1
## 23 ULAM 1
## 20 ROPS 3
## 6 BEPA 8
## 17 POGR4 9
## 18 PRSE2 9
## 1 ACRU 15
## 21 TIAM 16
## 11 OSVI 17
## 10 LARIX 26
## 19 QURU 34
## 16 PISY 37
## 4 BEAL2 51
## 7 FAGR 74
## 8 FRAM2 114
## 13 PIRE 131
## 12 PIAB 137
## 22 TSCA 253
## 15 PIST 277
## 2 ACSA3 349## Rename the columns
names(spptable) <- c("spp","freq")
barplot(spptable$freq,names.arg=spptable$spp,las=2)
### Exploring the aggregate function
## aggregate some of the data by plot
ba_plotsp <- aggregate(myData2$tree_ba,by=list(myData2$PLOT,myData2$tree_spp),sum)
head(ba_plotsp)## Group.1 Group.2 x
## 1 3 ACRU 1.069014
## 2 4 ACRU 1.310360
## 3 7 ACRU 2.821325
## 4 18 ACRU 1.115320
## 5 81 ACRU 1.993384
## 6 85 ACRU 1.082322## Reorder
ba_plotsp <- ba_plotsp[order(ba_plotsp$Group.1,ba_plotsp$Group.2),]
names(ba_plotsp) <- c("plot","spp","sba")## Aggregate stand totals
ba_plot <- aggregate(myData2$tree_ba,by=list(myData2$PLOT),sum)
names(ba_plot) <- c("plot","tba")
hist(ba_plot$tba,xlab="BA (m^2/acre)",nclass=40,col="blue")
head(ba_plot)## plot tba
## 1 1 22.847996
## 2 2 7.979427
## 3 3 20.434860
## 4 4 15.950673
## 5 5 8.527133
## 6 6 11.908818Subset to portion of data
####################################
# Subsett function on each stand
####################################
# define expansion factor
expf<-10
# Load Libraries
library(plyr)# STAND 1
stand1 <- subset(myData2, c(STAND==1))
#str(stand1)
table(stand1$tree_spp)##
## ACRU ACSA3 AMELA BEAL2 BELE BEPA FAGR FRAM2 JUCI LARIX OSVI PIAB
## 0 9 0 0 0 0 0 4 0 24 0 0
## PIRE PIRU PIST PISY POGR4 PRSE2 QURU ROPS TIAM TSCA ULAM
## 0 0 0 0 0 0 4 0 0 0 0boxplot(myData2$tree_dbh, xlab="Stand 1", ylab="DBH inches", col="blue")
summary(stand1$tree_dbh)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.5 12.3 16.6 15.5 18.7 24.1#--------------------# STAND 2
stand2 <- subset(myData2, c(STAND==2))
#str(stand2)
table(stand2$tree_spp)##
## ACRU ACSA3 AMELA BEAL2 BELE BEPA FAGR FRAM2 JUCI LARIX OSVI PIAB
## 0 0 0 0 0 0 0 0 0 0 0 1
## PIRE PIRU PIST PISY POGR4 PRSE2 QURU ROPS TIAM TSCA ULAM
## 0 0 23 0 0 0 0 0 0 0 0boxplot(stand2$tree_dbh, xlab="Stand 2", ylab="DBH inches", col="blue")
summary(stand2$tree_dbh)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.70 22.95 25.20 25.34 27.45 35.80summary(stand2$tree_ba) # per tree (ft^2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.6244 2.8730 3.4640 3.6240 4.1100 6.9900# Basal Area (ft^2)/acre
stand2<-subset(myData2, c(STAND==2))
#head(stand2)
stand2$ba_acre <-(length(stand2$STAND)/length(unique(stand2$PLOT)))*expf
stand2BA<-unique(stand2$ba_acre)
print(stand2BA)## [1] 120# Tree Per Acre
stand2$TPA<-sum(expf/stand2$tree_ba)/length(unique(stand2$PLOT))
stand2TPA <- unique(stand2$TPA)
stand2TPA## [1] 41.00533More with the plyr library
####################################################
## using the plyr library
#############################################
library(plyr)
#install.packages("Rmisc")
library(Rmisc)## Loading required package: lattice# Species richness by stand
sppR2<-ddply(myData2,"STAND",
summarize,
richness=length(unique(tree_spp)))
barplot(height=sppR2$richness, names.arg = sppR2$STAND,xlab="Stand #", ylab="Species Richness", col="blue")
### More to come