This runs you through the output of the 'best' model as it stands
rm(list = ls())
options(rpubs.upload.method = "internal")
options(warn = 0)
library(rgdal)
## Loading required package: sp
## rgdal: version: 0.8-8, (SVN revision 463) Geospatial Data Abstraction
## Library extensions to R successfully loaded Loaded GDAL runtime: GDAL
## 1.9.2, released 2012/10/08 Path to GDAL shared files:
## C:/Users/mmann/Documents/R/win-library/2.15/rgdal/gdal GDAL does not use
## iconv for recoding strings. Loaded PROJ.4 runtime: Rel. 4.7.1, 23
## September 2009, [PJ_VERSION: 470] Path to PROJ.4 shared files:
## C:/Users/mmann/Documents/R/win-library/2.15/rgdal/proj
library(car)
## Loading required package: MASS
## Loading required package: nnet
library(raster)
## Attaching package: 'raster'
## The following object(s) are masked from 'package:MASS':
##
## area, select
library(snow)
library(dismo)
library(pscl)
## Loading required package: mvtnorm
## Loading required package: coda
## Loading required package: lattice
## Loading required package: gam
## Loading required package: splines
## Loaded gam 1.06.2
## Loading required package: vcd
## Loading required package: grid
## Loading required package: colorspace
## Attaching package: 'vcd'
## The following object(s) are masked from 'package:raster':
##
## mosaic
## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis
library(maptools)
## Loading required package: foreign
## Checking rgeos availability: FALSE Note: when rgeos is not available,
## polygon geometry computations in maptools depend on gpclib, which has a
## restricted licence. It is disabled by default; to enable gpclib, type
## gpclibPermit()
library(fields)
## Loading required package: spam
## Spam version 0.29-2 (2012-08-17) is loaded. Type 'help( Spam)' or 'demo(
## spam)' for a short introduction and overview of this package. Help for
## individual functions is also obtained by adding the suffix '.spam' to the
## function name, e.g. 'help( chol.spam)'.
## Attaching package: 'spam'
## The following object(s) are masked from 'package:base':
##
## backsolve, forwardsolve
## Loading required package: maps
library(knitr)
source("G:\\Faculty\\Mann\\Share\\Scripts\\ten_sample.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\ten_samplev2.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\ten_samplev3.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\outofsamplehist2.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\outofsamplehist4.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\cfac.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\fire_plotter.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\py.rasterize2aggregate.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\py.rasterize.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\py.resample.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\py.EucDistance.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\py.ProjectRaster_management.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\nonlinear_test_zeroinfl2.R")
source("G:\\Faculty\\Mann\\Share\\Scripts\\multi_nonlinear_test_zeroinfl.R")
load("G:\\Faculty\\Mann\\Other\\Fire\\StackOutputs\\firedata_allfire3e.RData")
data1951_1975$PubLandDum = data1951_1975$PubLand == 0
data1976_2000$PubLandDum = data1976_2000$PubLand == 0
data1951_1975$allstationdist2 = pmin(data1951_1975$stationdist2, data1951_1975$airdist2)
data1976_2000$allstationdist2 = pmin(data1976_2000$stationdist2, data1976_2000$airdist2)
data1951_1975$den3x3a = data1951_1975$den3x3
data1951_1975$den3x3a[data1951_1975$den != 0] = 0 # put 0 for populated areas
data1976_2000$den3x3a = data1976_2000$den3x3
data1976_2000$den3x3a[data1976_2000$den != 0] = 0 # put 0 for populated areas
data1951_1975$den5x5a = data1951_1975$den5x5
data1951_1975$den5x5a[data1951_1975$den != 0] = 0 # put 0 for populated areas
data1976_2000$den5x5a = data1976_2000$den5x5
data1976_2000$den5x5a[data1976_2000$den != 0] = 0 # put 0 for populated areas
data1951_1975$den10x10a = data1951_1975$den10x10
data1951_1975$den10x10a[data1951_1975$den != 0] = 0 # put 0 for populated areas
data1976_2000$den10x10a = data1976_2000$den10x10a
data1976_2000$den10x10a[data1976_2000$den != 0] = 0 # put 0 for populated areas
zeroinflate_it <- function(equations_list, data_in_z, back_cast_data_z, distribution) {
split_eq = strsplit(equations_list[[1]], "=")
equation_name <<- trim_white_space(split_eq[[1]][1])
equation = as.formula(split_eq[[1]][2])
assign(equation_name, pscl::zeroinfl(equation, dist = distribution, data = data1976_2000[data1976_2000$train ==
T, ]))
AICER = AIC(get(equation_name))
# compare actual and predicted (actual-predicted for most graphs)
listed_out <<- ten_samplev3(equation_in = equation, graph_name = "e3f",
num_runs = 30, data_in = data_in_z, back_cast_data = back_cast_data_z,
distribution = distribution, pre_zip = pscl::zeroinfl(equation, dist = distribution,
data = data_in_z[data_in_z$train == T, ]), extraplots = F)
ten_sample_zip <<- listed_out[[1]]
vcov_list <<- list(listed_out[[2]])
zip_list <<- listed_out[[3]]
# assign( 'ten_sample_zip' ,ten_sample_zip, envir=.GlobalEnv) #something
# going wrong here not always assigning
return(list(ten_sample_zip, vcov_list))
}
equations_list = list("mikenewd4a = fire_poly ~cwdAave2+I(cwdAave2^2)+ cwdAsd2+log(aetAave2)+cwdAsd2:log(aetAave2) +Slope +PubLandDum | campdist2+ den+I(den^2) +Light5yrWIS2 +aetAave2+aetAave2:Light5yrWIS2 ")
#'#############################################'
#'Variable Definitions'
#'cwdAave2 = average annual mean cwd for 25 year period'
#'cwdAave2 = average annual standev of cwd for 25 year period'
#'aetAave2 = average annual mean aet for 25 year period'
#'PubLandDum = 1 for all public lands'
#'campdist2 = distance in meters for any campsite in CA listed in reserve america online listings'
#'den2 = housing density houses/ac'
#'I(den2^2) = (housing density houses/ac)^2'
#'Light5yrWIS = number of lightning strikes for most recent 5 year period provided by WIST'
#'aetAave2:Light5yrWIS2 = interaction term of lightning and aet'
#'#############################################'
listed_out2 = lapply(1:length(equations_list), function(x) zeroinflate_it(equations_list[[x]],
data_in_z = data1976_2000, back_cast_data_z = data1951_1975, distribution = "negbin"))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: NaNs produced
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: unimplemented pch value '26'
## Warning: unimplemented pch value '26'
## Warning: unimplemented pch value '26'
## Warning: unimplemented pch value '26'
## Warning: unimplemented pch value '26'
## Warning: unimplemented pch value '26'
## Warning: unimplemented pch value '26'
## Warning: unimplemented pch value '26'
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: unimplemented pch value '27'
## Warning: unimplemented pch value '27'
## Warning: unimplemented pch value '27'
## Warning: unimplemented pch value '27'
## Warning: unimplemented pch value '27'
## Warning: unimplemented pch value '27'
## Warning: unimplemented pch value '27'
## Warning: unimplemented pch value '27'
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: unimplemented pch value '28'
## Warning: unimplemented pch value '28'
## Warning: unimplemented pch value '28'
## Warning: unimplemented pch value '28'
## Warning: unimplemented pch value '28'
## Warning: unimplemented pch value '28'
## Warning: unimplemented pch value '28'
## Warning: unimplemented pch value '28'
## Warning: unimplemented pch value '29'
## Warning: unimplemented pch value '29'
## Warning: unimplemented pch value '29'
## Warning: unimplemented pch value '29'
## Warning: unimplemented pch value '29'
## Warning: unimplemented pch value '29'
## Warning: unimplemented pch value '29'
## Warning: unimplemented pch value '29'
## Warning: unimplemented pch value '30'
## Warning: unimplemented pch value '30'
## Warning: unimplemented pch value '30'
## Warning: unimplemented pch value '30'
## Warning: unimplemented pch value '30'
## Warning: unimplemented pch value '30'
## Warning: unimplemented pch value '30'
## Warning: unimplemented pch value '30'
## [1] "mean coefficients and pvalues"
## estimate p_val N_Sig
## (Intercept) -4.80238 1.791e-19 30
## cwdAave2 7.11398 2.191e-04 30
## I(cwdAave2^2) -2.74152 1.022e-02 29
## cwdAsd2 -1.92548 6.991e-02 26
## log(aetAave2) -1.04471 1.640e-01 19
## Slope 0.02483 5.899e-05 30
## PubLandDumTRUE 0.49123 1.166e-04 30
## cwdAsd2:log(aetAave2) 2.19853 8.628e-02 26
## Theta 1.79943 9.569e-02 26
## estimate p_val N_Sig
## (Intercept) -4.5961 0.106186 24
## campdist2 1.4883 0.002216 30
## den 3.5732 0.173401 21
## I(den^2) -0.2108 0.363616 10
## Light5yrWIS2 7.1067 0.037934 27
## aetAave2 0.2588 0.110687 22
## Light5yrWIS2:aetAave2 -8.3241 0.094326 23
## [1] "The mean results of the 15.5 regressions are as follows:"
multi_sample_zip = listed_out2[[1]][[1]]
vcov_list = listed_out2[[1]][[2]]
#'#############################################'
#'Plot Upper: % difference between predicted and actual counts red is mean model of 30 runs'
#'Plot lower: % difference between predicted and actual counts for backcast 1951-1975'
#'#############################################'
outofsamplehist4(equation_in = equations_list[[1]], zip = multi_sample_zip,
data_in_z = data1976_2000)
##
## Call:
## pscl::zeroinfl(formula = equation, data = data_in_z[data_in_z$train ==
## T, ], dist = distribution)
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.7892 -0.3344 -0.1986 -0.0858 22.4934
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.80238 0.67192 -7.15 8.9e-13 ***
## cwdAave2 7.11398 1.49631 4.75 2.0e-06 ***
## I(cwdAave2^2) -2.74152 0.77140 -3.55 0.00038 ***
## cwdAsd2 -1.92548 0.69297 -2.78 0.00546 **
## log(aetAave2) -1.04471 1.61185 -0.65 0.51689
## Slope 0.02483 0.00538 4.61 4.0e-06 ***
## PubLandDumTRUE 0.49123 0.09642 5.09 3.5e-07 ***
## cwdAsd2:log(aetAave2) 2.19853 1.58513 1.39 0.16545
## Log(theta) 1.36940 0.63997 2.14 0.03237 *
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.5961 2.0141 -2.28 0.02249 *
## campdist2 1.4883 0.3572 4.17 3.1e-05 ***
## den 3.5732 0.9675 3.69 0.00022 ***
## I(den^2) -0.2108 0.0429 -4.91 9.1e-07 ***
## Light5yrWIS2 7.1067 0.8583 8.28 < 2e-16 ***
## aetAave2 0.2588 1.1047 0.23 0.81476
## Light5yrWIS2:aetAave2 -8.3241 0.5783 -14.39 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 3.933
## Number of iterations in BFGS optimization: 49
## Log-likelihood: -1.93e+03 on 16 Df
## estobs0to3
## [1,] 271109.78
## [2,] 20113.31
## [3,] 1853.04
## [4,] 175.75
## [5,] 17.91
#'#############################################'
#'Plot upper left: Probability of 1 or more fires'
#'Plot upper right: mean fire return interval'
#'Plot lower left: Expected number of fires'
#'Plot lower right: bar = actual # pixels with 1-5 fire events, triangle = predicted number of pixels with 1-5 fire events'
#'#############################################'
#'#############################################'
#'Joint F tests for nonlinear variables and interactions'
#'#############################################'
multi_nonlinear_test_zeroinfl(variable_name_in = "cwdAave2", root_in = 2, zip_list_in = zip_list,
count_or_zero_in = "count")
## F Pr(>F)
## 1 34.40 4.674e-09
## 2 62.59 2.915e-15
## 3 50.13 1.575e-12
## 4 52.61 4.494e-13
## 5 64.21 1.295e-15
## 6 29.07 7.205e-08
## 7 53.72 2.556e-13
## 8 37.20 1.118e-09
## 9 26.02 3.472e-07
## 10 59.55 1.351e-14
## 11 51.92 6.350e-13
## 12 47.90 4.873e-12
## 13 59.74 1.226e-14
## 14 56.90 5.144e-14
## 15 27.21 1.874e-07
## 16 31.50 2.069e-08
## 17 43.44 4.669e-11
## 18 51.49 7.898e-13
## 19 21.92 2.887e-06
## 20 59.82 1.177e-14
## 21 37.37 1.027e-09
## 22 33.22 8.551e-09
## 23 35.82 2.272e-09
## 24 21.75 3.157e-06
## 25 63.38 1.965e-15
## 26 57.50 3.791e-14
## 27 29.26 6.540e-08
## 28 55.42 1.083e-13
## 29 56.60 5.974e-14
## 30 27.56 1.563e-07
## [1] "# significant at 5%: 30"
## [1] "# significant at 10%: 30"
## [1] "mean F value: 44.6410643661274"
## [1] "mean F P-value: 2.30362432703625e-07"
multi_nonlinear_test_zeroinfl(variable_name_in = "den", root_in = 2, zip_list_in = zip_list,
count_or_zero_in = "zero")
## F Pr(>F)
## 1 14.882052 1.154e-04
## 2 10.486978 1.208e-03
## 3 1.806151 1.790e-01
## 4 16.701616 4.421e-05
## 5 1.076794 2.995e-01
## 6 5.559460 1.841e-02
## 7 3.556312 5.936e-02
## 8 18.023015 2.209e-05
## 9 0.441175 5.066e-01
## 10 15.735769 7.353e-05
## 11 0.169418 6.806e-01
## 12 0.111166 7.388e-01
## 13 1.578065 2.091e-01
## 14 8.822939 2.984e-03
## 15 15.481515 8.409e-05
## 16 3.956009 4.674e-02
## 17 1.392387 2.380e-01
## 18 19.840511 8.543e-06
## 19 6.764448 9.318e-03
## 20 16.292800 5.482e-05
## 21 5.513878 1.889e-02
## 22 8.415532 3.731e-03
## 23 4.620385 3.163e-02
## 24 0.142763 7.056e-01
## 25 19.063297 1.282e-05
## 26 14.934226 1.123e-04
## 27 12.191269 4.830e-04
## 28 0.002523 9.599e-01
## 29 15.949012 6.571e-05
## 30 15.125457 1.015e-04
## [1] "# significant at 5%: 20"
## [1] "# significant at 10%: 21"
## [1] "mean F value: 8.62123069091938"
## [1] "mean F P-value: 0.157019301503275"
multi_nonlinear_test_zeroinfl(variable_name_in = c("cwdAsd2", "log(aetAave2)",
"cwdAsd2:log(aetAave2)"), root_in = 0, zip_list_in = zip_list, count_or_zero_in = "count")
## F Pr(>F)
## 1 47.226 3.201e-30
## 2 18.909 3.264e-12
## 3 9.775 1.967e-06
## 4 15.430 5.276e-10
## 5 9.420 3.285e-06
## 6 4.102 6.430e-03
## 7 9.870 1.714e-06
## 8 3.393 1.714e-02
## 9 11.830 1.000e-07
## 10 7.258 7.379e-05
## 11 1.758 1.529e-01
## 12 11.964 8.230e-08
## 13 16.105 1.969e-10
## 14 14.382 2.434e-09
## 15 7.512 5.124e-05
## 16 11.875 9.365e-08
## 17 12.151 6.265e-08
## 18 5.974 4.612e-04
## 19 3.734 1.072e-02
## 20 10.387 8.112e-07
## 21 14.870 1.194e-09
## 22 4.197 5.639e-03
## 23 13.873 5.111e-09
## 24 9.526 2.818e-06
## 25 11.466 1.697e-07
## 26 19.567 1.248e-12
## 27 10.231 1.016e-06
## 28 5.232 1.319e-03
## 29 9.561 2.679e-06
## 30 8.890 7.058e-06
## [1] "# significant at 5%: 29"
## [1] "# significant at 10%: 29"
## [1] "mean F value: 11.3489897708004"
## [1] "mean F P-value: 0.00649199795133487"
multi_nonlinear_test_zeroinfl(variable_name_in = c("Light5yrWIS2", "aetAave2",
"Light5yrWIS2:aetAave2"), root_in = 0, zip_list_in = zip_list, count_or_zero_in = "zero")
## F Pr(>F)
## 1 15.014 9.686e-10
## 2 5.869 5.348e-04
## 3 19.041 2.694e-12
## 4 18.829 3.669e-12
## 5 6.650 1.760e-04
## 6 12.930 2.017e-08
## 7 16.407 1.267e-10
## 8 17.011 5.240e-11
## 9 12.366 4.587e-08
## 10 21.366 8.948e-14
## 11 22.952 8.767e-15
## 12 16.088 2.018e-10
## 13 3.132 2.453e-02
## 14 12.580 3.362e-08
## 15 22.508 1.681e-14
## 16 19.401 1.590e-12
## 17 20.288 4.341e-13
## 18 18.932 3.158e-12
## 19 17.092 4.653e-11
## 20 14.746 1.431e-09
## 21 16.225 1.654e-10
## 22 15.187 7.527e-10
## 23 15.660 3.772e-10
## 24 14.828 1.271e-09
## 25 16.921 5.980e-11
## 26 14.816 1.292e-09
## 27 16.411 1.260e-10
## 28 18.681 4.557e-12
## 29 22.296 2.294e-14
## 30 13.829 5.454e-09
## [1] "# significant at 5%: 30"
## [1] "# significant at 10%: 30"
## [1] "mean F value: 15.9352176272228"
## [1] "mean F P-value: 0.000841486608413654"
############################################################## backcast
############################################################## plots hold
############################################################## natural and
############################################################## human
############################################################## variables
############################################################## at mean
############################################################## choose ##
backcast_period = "1951_1975"
estimation_period = "1976_2000" # period used in ZIP model above
sample_zip = multi_sample_zip # choose which estimated regression to use. to use second put 2
## run ## get data
backcast_data_z = get(paste("data", backcast_period, sep = ""), envir = .GlobalEnv)
current_data_z = get(paste("data", estimation_period, sep = ""), envir = .GlobalEnv)
# add extra variables created
backcast_data_z$allstationdist2 = pmin(data1951_1975$stationdist2, data1951_1975$airdist2)
current_data_z$allstationdist2 = pmin(data1976_2000$stationdist2, data1976_2000$airdist2)
# predict using mean coefficients from zeroinflate_it backcast fire
# probability of no fire
backcast = predict(sample_zip, backcast_data_z, type = "prob")
backcast = backcast[, 1] # prob of zero fires
raster_example = raster(universal_variables[1]) # example raster to use in rasterize
pts = backcast_data_z[, c("x", "y")] # points that match backcast_data
ptscurrent = current_data_z[, c("x", "y")] # points that match backcast_data
# then
backcast_prob_no_fire <- rasterize(pts, raster_example, backcast, fun = sum)
backcast_prob_atleastone_fire <<- 1 - backcast_prob_no_fire
windows()
plot(backcast_prob_atleastone_fire, main = "probability of 1 fire backcast")
# current fire probability of no fire
current_prob_fire = predict(sample_zip, current_data_z, type = "prob")
current_prob_fire = current_prob_fire[, 1]
current_prob_no_fire <- rasterize(ptscurrent, raster_example, current_prob_fire,
fun = sum)
current_prob_atleastone_fire <<- 1 - current_prob_no_fire
windows()
plot(current_prob_atleastone_fire, main = "current prob at least one fire")
writeRaster(current_prob_atleastone_fire, filename = "C:/Users/mmann1123/Desktop/Share/Fire/Fire_All_Script_Outputs\\fire_prob_7599.tif",
format = "GTiff", overwrite = TRUE)
## Error: Attempting to write a file to a path that does not exist:
## C:/Users/mmann1123/Desktop/Share/Fire/Fire_All_Script_Outputs
# current fire count
current_resp_fire <- predict(sample_zip, current_data_z, type = "response")
current_resp_fire <<- rasterize(ptscurrent, raster_example, current_resp_fire,
fun = sum) # raster of human influence
writeRaster(current_resp_fire, filename = "C:/Users/mmann1123/Desktop/Share/Fire/Fire_All_Script_Outputs\\fire_count_7599.tif",
format = "GTiff", overwrite = TRUE)
## Error: Attempting to write a file to a path that does not exist:
## C:/Users/mmann1123/Desktop/Share/Fire/Fire_All_Script_Outputs
windows()
plot(current_resp_fire, main = "current expected fire count")
# backcast fire count
backcast_resp_fire <- predict(sample_zip, backcast_data_z, type = "response")
backcast_resp_fire <<- rasterize(pts, raster_example, backcast_resp_fire, fun = sum) # raster of human influence
windows()
plot(backcast_resp_fire, main = "Backcast expected fire count")
# difference in coutn
windows()
plot(current_resp_fire - backcast_resp_fire, main = "Difference in count back and now")
# mean FRI
windows()
library(fields)
image.plot(25/current_resp_fire, main = "current Mean Fire return", zlim = c(0,
250))
mri = 25/current_resp_fire
# writeRaster(mri,
# filename='C:/Users/mmann1123/Desktop/Share/Fire/Fire_All_Script_Outputs\\fire_FRI_7599.tif',
# format='GTiff', overwrite=TRUE)
human_var_names_to_search = c("den", "den3x3", "den5x5", "den10x10", "NPark",
"PubLand", "AllRoads", "Incorp", "PrimSec", "Pp30k", "Pp20k", "Ppall", "campdist",
"stationdist", "airdist") # not using grep, be explicit
public_land = c("PubLandDum")
human_wo_public_var_names_to_search = human_var_names_to_search %w/o% public_land
factors = c()
sample_zip = multi_sample_zip
# get mean values for current
current_mn = as.data.frame(t(colMeans(current_data_z, na.rm = T)))
current_names = names(current_mn)
rowsc = dim(current_data_z)[1]
current_mn = matrix(rep(as.numeric(current_mn), rowsc), rowsc, byrow = T) # createdataframe of all column means while matching dim of input data
colnames(current_mn) = current_names
human_influence_loc_current = multi_grep_character(human_var_names_to_search,
colnames(current_mn)) #col # of human variables
if (!is.null(factors)) {
human_influence_loc_current = unique(human_influence_loc_current %w/o% multi_grep_character(factors,
colnames(current_mn)))
} # remove factors so doesn't hold at mean
wo_human_influence_loc_current = seq(1, dim(current_mn)[2]) %w/o% human_influence_loc_current
if (!is.null(factors)) {
wo_human_influence_loc_current = unique(wo_human_influence_loc_current %w/o%
multi_grep_character(factors, colnames(current_mn)))
} # remove factors so doesn't hold at mean
null = as.data.frame(current_mn) # has to be data.frame here to set logicals
null[, paste(public_land)] = F # set public land to 'false'
# hold all but natural at mean public land at zero
hm_mn_pb_0 = current_data_z # actual data
hm_mn_pb_0[, unique(human_influence_loc_current)] = current_mn[, unique(human_influence_loc_current)] # set human effects at mean values
hm_mn_pb_0[, paste(public_land)] = F # set public land 'off'
# turn 'on' public lands
hm_mn_pb_1 = hm_mn_pb_0
for (hold in public_land) {
# use this, don't switch all to 1, b/c that would be as if all lands were public
hm_mn_pb_1[, paste(hold)] = current_data_z[, paste(hold)]
}
# turn on all, public land =1 where public land exists
all_normal_pb_1 = current_data_z
null = predict(object = sample_zip, newdata = as.data.frame(null), type = "response") # predict for all held at mean and publand=0
null = rasterize(ptscurrent, raster_example, null, fun = sum) # convert to matrix
hm_mn_pb_0 = predict(sample_zip, hm_mn_pb_0, type = "response")
hm_mn_pb_0 = rasterize(ptscurrent, raster_example, hm_mn_pb_0, fun = sum)
hm_mn_pb_1 = predict(sample_zip, hm_mn_pb_1, type = "response")
hm_mn_pb_1 = rasterize(ptscurrent, raster_example, hm_mn_pb_1, fun = sum)
all_normal_pb_1 = predict(sample_zip, all_normal_pb_1, type = "response")
all_normal_pb_1 = rasterize(ptscurrent, raster_example, all_normal_pb_1, fun = sum)
pub_effect = (hm_mn_pb_1 - hm_mn_pb_0)
nat_effect = (hm_mn_pb_0) # - null
hum_effect = (all_normal_pb_1 - hm_mn_pb_1)
color = terrain.colors(45)
color1 = color
color1[1] = "#BDBDBD"
plot((nat_effect), main = "natural effects", col = color1)
plot((pub_effect), main = " Public land effects", col = color1, zlim = c(0,
maxValue(nat_effect)))
color2 = color[length(color):1]
plot(hum_effect, main = " Human effects", col = color2)
