Something about how I've written this optimize function is not letting it refer to existing values in my "Input" data set to fill in a new column, "Opt":
Input=read.csv("....csv")
Input$Opt=0
Input$Opt <- optimize(f = function(x) abs(10.16 - KozakTaper(Bark='ob',
SPP=Input$Species,
DHT=x,
DBH=Input$DBH,
HT=Input$Ht,
Planted=0)),
lower=Input$Ht*.25, upper=Input$Ht+1,
maximum = FALSE, tol = .Machine$double.eps^0.25)[[1]]
I get the error
invalid function value in 'optimize'
Here is a brief definition of "KozakTaper" so you have some context for what I'm trying to do. The trick, and why I need to use optimize or something like it, is that I don't want y (which is what KozakTaper returns). I want to know what DHT is if y=10.16, but because the equation can't be rearranged to solve for DHT, I'm using optimize to return the DHT value that minimizes the difference between y and 10.16.
KozakTaper=function(Bark,SPP,DHT,DBH,HT,Planted){
if(Bark=='ob' & SPP=='AB'){
a0_tap=1.0693567631
a1_tap=0.9975021951
a2_tap=-0.01282775
b1_tap=0.3921013594
b2_tap=-1.054622304
b3_tap=0.7758393514
b4_tap=4.1034897617
b5_tap=0.1185960455
b6_tap=-1.080697381
b7_tap=0}
else if(Bark=='ob' & SPP=='RS'){
a0_tap=0.8758
a1_tap=0.992
a2_tap=0.0633
b1_tap=0.4128
b2_tap=-0.6877
b3_tap=0.4413
b4_tap=1.1818
b5_tap=0.1131
b6_tap=-0.4356
b7_tap=0.1042}
else{
a0_tap=1.1263776728
a1_tap=0.9485083275
a2_tap=0.0371321602
b1_tap=0.7662525552
b2_tap=-0.028147685
b3_tap=0.2334044323
b4_tap=4.8569609081
b5_tap=0.0753180483
b6_tap=-0.205052535
b7_tap=0}
p = 1.3/HT
z = DHT/HT
Xi = (1 - z^(1/3))/(1 - p^(1/3))
Qi = 1 - z^(1/3)
y = (a0_tap * (DBH^a1_tap) * (HT^a2_tap)) * Xi^(b1_tap * z^4 + b2_tap * (exp(-DBH/HT)) +
b3_tap * Xi^0.1 + b4_tap * (1/DBH) + b5_tap * HT^Qi + b6_tap * Xi + b7_tap*Planted)
return(y=round(y,4))}
If you have a different approach for finding DHT for each row of data, I'm open to other suggestions. The optimize strategy works fine for one data point (a tree of given species, diameter, and height), but since I've put in the references to column names, it's hitting the error. All advice on what to modify is appreciated!
-KB
Research Fellow, Appalachian Mountain Club
A simplified example of the input data:
> dput(head(Input))
structure(list(Species = structure(c(3L, 3L, 3L, 3L, 8L, 8L), .Label = c("AB",
"BC", "BF", "EH", "PB", "PR", "RM", "RS", "SM", "ST", "WA", "WP",
"YB"), class = "factor"), DBH = c(6.9000001, 8.1000004, 5.8000002,
6.5999999, 9.5, 7.5999999), Ht = c(44, 43, 34, 41, 56, 58)), .Names = c("Species",
"DBH", "Ht"), row.names = c(NA, 6L), class = "data.frame")