##compute and graph the frequency histogram of tree population density #m is the number of column of census in data1 # for 5 m x 5 m scale newfielddata=read.csv("newfielddata.csv",header=T) sitegroup=read.csv("sitegroup.csv",header=T) site=read.csv("site.csv",header=T) Moments1=function(data1,m,censusyear="1974") { library(PearsonDS) a1=na.omit(data1[,c(2,3,m)]) #draw the quadrat,sbuplot and dbh informaton for focus cencus a1=a1[a1[,3]>=30,] # remove DBH<30mm trees # for 5 m x 5 m scale density25=vector() #generate a empty vector of mean density (no/ha) of each 5m x 5 m subplot density50.1=vector()#generate a empty vector of mean density (no/ha) of each 5m x 10 m subplot density50.2=vector() #generate a empty vector of mean density (no/ha) of each 10m x 5 m subplot density400=vector() #generate a empty vector of mean density (no/ha) of each 20m x 20 m subplot qrat=sort(unique(a1$quadrat)) # extract unique quadrat, then order the vector qn=length(qrat) # compute number of unique quadrat quadrat=rep(qrat,rep(16,qn)) # generate a new vector for unique quadrat, each quadrat has 16 subplot subplot=rep(1:16,qn) # generate a new vector for subplot, each subplot was repeated qn times n=qn*16 # compute the lenght of variance and mean biomass, equal to total count of all 5 m x 5 m subplots in data Density=vector() # generate a empty vector of mean density (no/ha) of each 5 m x 5 m subplot for(i in 1:n) { Density[i]=length(a1[a1$quadrat==quadrat[i]&a1$subplot==subplot[i],]$subplot)*400 } density25=Density #compute expected slope according to the above formulas of Cohen and Xu (unpublished) em1=empMoments(density25) expslope1=em1[[1]]*em1[[3]]/(em1[[2]])^0.5 n1=length(density25) N=n1 Mean=em1[[1]] Variance=em1[[2]] Skewness=em1[[3]] Kurtosis=em1[[4]] hist((density25),main="", freq=T,xlab="Density (#/ha)",ylab="Frequency") title(main="5 m x 5 m",line=1) title(main=censusyear,line=2.5) legend("topright",c(paste("N:", round(n1,0)),paste("Mean:", round(em1[[1]],0)), paste("Variance:", round(em1[[2]],0)),paste("Skewness:", round(em1[[3]],2)), paste("Kurtosis:", round(em1[[4]],2))),bty="n",cex = 1.3) result=data.frame(N,Mean,Variance,Skewness,Kurtosis) return(result) } par(mfrow=c(2,5)) Moments1(newfielddata,4) Moments(newfielddata,5,censusyear="1976") Moments1(newfielddata,6,censusyear="1978") Moments1(newfielddata,7,censusyear="1980") Moments1(newfielddata,8,censusyear="1983") Moments1(newfielddata,9,censusyear="1985") Moments1(newfielddata,10,censusyear="1987") Moments1(newfielddata,11,censusyear="1993") Moments1(newfielddata,12,censusyear="1999") Moments1(newfielddata,13,censusyear="2004") ##compute and graph the frequency histogram of tree population density #m is the number of column of census in data1 # for 20 m x 20 m scale Moments2=function(data1,m,censusyear="1974") { library(PearsonDS) a1=na.omit(data1[,c(2,3,m)]) #draw the quadrat,sbuplot and dbh informaton for focus cencus a1=a1[a1[,3]>=30,] # remove DBH<30mm trees # for 20 m x 20 m scale density400=table(a1$quadrat)*25 #compute expected slope according to the above formulas of Cohen and Xu (unpublished) em4=empMoments(density400) expslope4=em4[[1]]*em4[[3]]/(em4[[2]])^0.5 n4=length(density400) N=n4 Mean=em4[[1]] Variance=em4[[2]] Skewness=em4[[3]] Kurtosis=em4[[4]] hist((density400),main="", freq=T,xlab="Density (#/ha)",ylab="Frequency") title(main="20 m x 20 m",line=1) title(main=censusyear,line=2.5) legend("topright",c(paste("N:", round(n4,0)),paste("Mean:", round(em4[[1]],0)), paste("Variance:", round(em4[[2]],0)),paste("Skewness:", round(em4[[3]],2)), paste("Kurtosis:", round(em4[[4]],2))),bty="n",cex = 1.3) result=data.frame(N,Mean,Variance,Skewness,Kurtosis) return(result) } Moments2(newfielddata,4) Moments2(newfielddata,5,censusyear="1976") Moments2(newfielddata,6,censusyear="1978") Moments2(newfielddata,7,censusyear="1980") Moments2(newfielddata,8,censusyear="1983") Moments2(newfielddata,9,censusyear="1985") Moments2(newfielddata,10,censusyear="1987") Moments2(newfielddata,11,censusyear="1993") Moments2(newfielddata,12,censusyear="1999") Moments2(newfielddata,13,censusyear="2004")