Nevins_Homwork_9
Matthias Nevins
March 28, 2017
myData<-read.csv("MBR_Overstory_obs_wrking_demo.csv", header=T)
#str(myData)
dbh <- myData$tree_dbh
ba <- myData$tree_ba
stand <- myData$STAND
species <- myData$tree_spp
data <- data.frame(stand, species,dbh,ba)
str(data)## 'data.frame': 1666 obs. of 4 variables:
## $ stand : int 1 1 1 1 1 1 1 1 1 1 ...
## $ species: Factor w/ 24 levels "ACRU","ACSA3",..: 10 2 2 10 2 2 10 10 2 10 ...
## $ dbh : num 18.7 12.3 8.8 24.1 13.7 7.1 15.4 15.6 8.2 13.9 ...
## $ ba : num 1.907 0.825 0.422 3.168 1.024 ...#print(data)REGRESSION ANALYSIS
###########################################################################
# FUNCTION: Regression
# input: x = 1:10, y = runif(10)
# output: linear model list and data summary
# This function runs a basic linear regression
####################################################################
Regression <- function(x=1:10, y=runif(10)){
rmod <- lm(y~x)
out <- c(slope=summary(rmod)$coefficients[2,1],
pValue=summary(rmod)$coefficients[2,4],
tvalue=summary(rmod)$coefficients[2,3],
standarderror=summary(rmod)$coefficients[2,2],
RSquared=summary(rmod)$adj.r.squared)
return(out)
}
#####################################################################
# END OF FUNCITON
#####################################################################Test Function
Regression()## slope pValue tvalue standarderror RSquared
## 0.04176821 0.26309817 1.20373575 0.03469882 0.04751622myData Run
Regression(x=data$dbh, y=data$ba)## slope pValue tvalue standarderror RSquared
## 2.370699e-01 0.000000e+00 1.315100e+02 1.802676e-03 9.121783e-01###########################################################################
# FUNCTION: lrPlot
# lrPlot() plots a basic linear regression from the Regression() function
# input: x = 1:10, y = runif(10)
# output: linear regression plot
###########################################################################
lrPlot <- function(x=1:10, y=runif(10)){
rmod <- lm(y~x)
plot <- c(plot(x=x, y=y,pch=21, bg="lightblue", cex=1.5, main="LINEAR REGRESSION PLOT", xlab = "Tree DBH (inches)", ylab = "Tree Basal Area (sq. ft./acre"), abline(rmod))
return(plot)
}
###########################################################################
# END OF FUNCITON
###########################################################################TEST lrPlot()
lrPlot()
## NULLTest lrPlot on myData
lrPlot(x=data$dbh, y=data$ba)
## NULLANOVA
###########################################################################
# FUNCTION: ANOVA
# input: x = discrete variables (x= Acer, Pinus); y=runif(10))
# output: anova model list and data summary
####################################################################
ANOVA <- function(x=(c(rep("Acer", 3), rep("Pine", 7))), y=runif(10)){
amod <- aov(y~x)
out <- list(amod, summary(amod))
return(out)
}
#####################################################################
# END OF FUNCITON
#####################################################################Test Function
ANOVA()## [[1]]
## Call:
## aov(formula = y ~ x)
##
## Terms:
## x Residuals
## Sum of Squares 0.0713357 0.5094573
## Deg. of Freedom 1 8
##
## Residual standard error: 0.2523532
## Estimated effects may be unbalanced
##
## [[2]]
## Df Sum Sq Mean Sq F value Pr(>F)
## x 1 0.0713 0.07134 1.12 0.321
## Residuals 8 0.5095 0.06368myData Run
ANOVA(x=data$species, y=data$dbh)## [[1]]
## Call:
## aov(formula = y ~ x)
##
## Terms:
## x Residuals
## Sum of Squares 22694.86 81542.42
## Deg. of Freedom 23 1642
##
## Residual standard error: 7.047015
## Estimated effects may be unbalanced
##
## [[2]]
## Df Sum Sq Mean Sq F value Pr(>F)
## x 23 22695 986.7 19.87 <2e-16 ***
## Residuals 1642 81542 49.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ANOVA PLOT
###########################################################################
# FUNCTION: anovaPlot
# anovaPlot() plots anova analysis
# input: x=as.factor(c(rep("Acer", 3), rep("Pine", 7)), y=runif(10)
# output: anova plot
###########################################################################
anovaPLot <- function(x=as.factor(c(rep("Acer", 3), rep("Pine", 7))), y=runif(10)){
plot <- boxplot(y~x, main="ANOVA PLOT", xlab="Stand", ylab="Tree DBH - inches")
return(plot)
}
###########################################################################
# END OF FUNCITON
###########################################################################Test anovaPlot()
anovaPLot()
## $stats
## [,1] [,2]
## [1,] 0.03747161 0.2552529
## [2,] 0.04953616 0.3723352
## [3,] 0.06160071 0.5616851
## [4,] 0.46793552 0.5857273
## [5,] 0.87427034 0.7148026
##
## $n
## [1] 3 7
##
## $conf
## [,1] [,2]
## [1,] -0.3200688 0.4342508
## [2,] 0.4432702 0.6891194
##
## $out
## numeric(0)
##
## $group
## numeric(0)
##
## $names
## [1] "Acer" "Pine"anovaPlot(mydata)
anovaPLot(x=data$stand, y=data$dbh)
## $stats
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 6.5 19.5 8.10 9.10 5.20 11.3 4.00 6.00 5.1 5.5 5.40 11.5
## [2,] 12.0 21.3 11.10 11.05 8.00 14.8 12.35 8.95 7.5 7.2 8.45 11.5
## [3,] 16.5 24.1 16.15 12.25 15.05 17.6 15.30 13.50 9.6 12.1 15.60 14.8
## [4,] 18.7 27.0 18.60 13.35 23.60 23.4 18.40 19.55 11.5 34.2 22.20 15.1
## [5,] 24.1 33.6 26.70 14.70 30.40 33.5 25.80 24.60 11.6 46.0 34.90 16.6
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
## [1,] 9.20 4.20 15.6 14.50 12.30 13.70 7.70 8.10 4.8 12.80 8.20
## [2,] 13.15 8.05 17.3 15.15 14.85 20.85 12.15 10.55 12.4 15.10 15.05
## [3,] 15.05 14.35 19.5 16.50 16.00 23.50 13.65 13.10 15.2 16.45 17.30
## [4,] 18.00 16.25 23.1 16.85 17.45 26.15 15.85 17.65 18.0 18.80 20.20
## [5,] 22.50 28.40 25.5 18.10 20.50 30.80 19.80 22.10 24.7 22.20 22.00
## [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34]
## [1,] 8.60 5.60 14.2 6.60 14.50 5.20 5.1 5.50 4.50 4.60 6.20
## [2,] 13.25 8.55 15.3 12.10 16.60 11.70 8.3 10.55 12.15 10.40 17.35
## [3,] 15.30 14.00 16.4 16.55 18.45 16.10 10.4 18.90 15.30 14.70 22.20
## [4,] 18.15 16.60 17.4 26.20 19.95 22.75 14.8 23.75 18.75 23.85 24.80
## [5,] 25.10 26.90 19.3 32.60 22.10 32.90 24.4 29.00 24.90 38.90 28.40
## [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## [1,] 4.10 5.60 3.2 6.90 13.10 11.30 10.2 5.00 7.8 6.10
## [2,] 11.55 14.25 10.7 12.20 19.55 17.30 19.2 10.40 17.3 17.60
## [3,] 18.20 18.75 15.2 16.85 25.90 20.80 21.2 14.25 22.4 24.10
## [4,] 22.70 20.75 21.6 23.20 29.50 24.55 27.4 19.80 31.2 31.55
## [5,] 36.40 26.30 36.7 33.10 42.90 34.30 36.7 33.00 48.5 48.40
##
## $n
## [1] 43 33 30 60 10 35 47 35 16 11 20 5 40 16 17 11 63
## [18] 68 16 31 73 20 15 35 27 33 16 16 43 61 20 52 56 19
## [35] 43 48 194 42 15 35 15 50 56 75
##
## $conf
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 14.88565 22.53226 13.9865 11.78085 7.255618 15.30321 13.90568
## [2,] 18.11435 25.66774 18.3135 12.71915 22.844382 19.89679 16.69432
## [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
## [1,] 10.66907 8.02 -0.762474 10.74214 12.25625 13.83837 11.111 17.2774
## [2,] 16.33093 11.18 24.962474 20.45786 17.34375 16.26163 17.589 21.7226
## [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
## [1,] 15.69014 15.48244 22.4845 12.1885 11.08519 14.16442 15.14279 15.19904
## [2,] 17.30986 16.51756 24.5155 15.1115 15.11481 16.23558 17.75721 19.40096
## [,24] [,25] [,26] [,27] [,28] [,29] [,30]
## [1,] 13.99136 11.55223 15.82241 10.9805 17.12675 13.43753 9.085061
## [2,] 16.60864 16.44777 16.97759 22.1195 19.77325 18.76247 11.714939
## [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## [1,] 14.23646 13.8539 11.86022 19.49955 15.51343 17.26765 13.96353
## [2,] 23.56354 16.7461 17.53978 24.90045 20.88657 20.23235 16.43647
## [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## [1,] 14.16821 21.84086 18.86375 17.85478 12.14961 19.4652 21.55492
## [2,] 19.53179 29.95914 22.73625 24.54522 16.35039 25.3348 26.64508
##
## $out
## [1] 10.6 35.8 10.9 7.0 10.7 51.4 4.1 5.5 42.7 47.4 4.9 4.7 38.2 7.2
## [15] 47.6 4.7 32.3 34.5 56.1 52.0 34.7 52.3 48.8 44.0 35.8 61.4 49.4 18.8
## [29] 20.4 23.0 5.5 29.5 7.6 7.3 6.7 20.7 5.3 5.3 8.1 7.6 5.3 6.3
## [43] 11.3 34.4 8.8 8.8 34.5 12.8 4.3 4.0 31.1 35.8 26.5 27.0 25.3 4.3
## [57] 5.4 6.1 5.9 4.8 30.5 6.5 6.8 9.0 6.8 6.8 26.9 40.6 30.1 30.1
## [71] 38.4 30.7 31.6 32.2 31.1 36.5 34.2
##
## $group
## [1] 2 2 2 2 2 3 4 4 4 4 4 4 4 4 4 4 4 4 6 7 7 8 8
## [24] 8 8 8 8 9 9 9 12 13 15 15 16 16 17 17 17 17 17 18 18 18 18 18
## [47] 18 18 18 18 20 20 21 21 22 22 22 23 24 24 24 26 26 26 26 26 28 29 30
## [70] 30 32 32 32 32 36 42 42
##
## $names
## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14"
## [15] "15" "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28"
## [29] "30" "31" "32" "33" "34" "35" "37" "38" "39" "40" "41" "42" "43" "44"
## [43] "45" "46"library(fGarch)## Loading required package: timeDate## Loading required package: timeSeries## Loading required package: fBasics## ## Rmetrics Package fBasics## Analysing Markets and calculating Basic Statistics## Copyright (C) 2005-2014 Rmetrics Association Zurich## Educational Software for Financial Engineering and Computational Science## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.## https://www.rmetrics.org --- Mail to: info@rmetrics.orgMOCK DATA
# Assess true dbh
hist(dbh)
mean(dbh)## [1] 17.52137dbh2<-na.omit(dbh)
mean(dbh2)## [1] 17.52137sd(dbh2)## [1] 7.91233# mock data simulation
dbhNew<-rsnorm(1667, mean = 17.52, sd = 7.912, xi = 2.4)
hist(dbhNew)
# Assess true ba
hist(ba)
ba2<-na.omit(ba)
mean(ba2) #2.015## [1] 2.015668sd(ba2) #1.96## [1] 1.963941summary(ba2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.05585 0.78540 1.46700 2.01600 2.56800 20.56000#Evaluate rgamma distribution
hist(rgamma(n=100, shape = 1,scale = 10))
baNew <- rgamma(n=1667, shape = .2, scale = 10)
hist(baNew)
mean(baNew)## [1] 1.806252summary(baNew)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00596 0.19830 1.80600 1.53200 47.27000TEST REGRESSION FUNCTION ON MOCK DATA
# call to regression function in global envi
Regression(x=dbhNew, y=baNew)## slope pValue tvalue standarderror RSquared
## -0.003235539 0.802845014 -0.249706655 0.012957362 -0.000563130TEST lrPlot
# call to lrPlot function in global envi
lrPlot(x=dbhNew, y=baNew)
## NULL