Monte-Carlo Methods

Introduction

Monte-Carlo methods refer to running portions or all of the modeling process over and over again using random variation to help understand the nature and robustness of our models.

Monte-Carlo is often reduce to "MC" and sometimes, incorrlectly, to MCMC. MCMC stands for Monte-Carlo Marcov Chains and since we do not use Marcov Chains in the models on this web site, we'll just use MC.

Check out our section on loops if you are not familiary with loops in R.

MC in R

A simple example would be to take one of our GAM models and try different values for the gamma parameter. The code below will run through a range of gamma values and record AIC, Deviance, and the Effective Degrees of Freedom (EDF) for each model. Then, the percent deviance is computed and the range of the performance metrics are plotted for the range of gamma values. Run this code for one of your models and then select the "best" gamma vaule based on the results.

require(mgcv) # load the GAMs library
	
# load the file with the field measurements for DougFir height
TheData = read.csv("C:\\Temp\\DougFir_All_CA_8km.csv"

# Remove any NA (null) values
TheData=na.omit(TheData)

AICs=vector() # setup a vector to contain the AIC values
Deviances=vector() # setup a vector to contain the Deviance values
NullDeviances=vector() # setup a vector to contain the NullDeviance values
EDFs=vector() # setup a vector to contain the effective degrees of freedom values

# loop through the gamma values
Gammas=seq(1,10,by=1) # setup a sequence of values for gamma

Count=1
while (Count<=length(Gammas)) # go through all the gamma values
{
	# create a GAM model using a gamma value from our sequence
	TheGamma=Gammas[Count] # get the target gamma value from our sequence
	
	TheModel=gam(Height~s(MaxTemp)+s(MinTemp)+s(AnnualPrecip),data=TheData,gamma=TheGamma)
	
	# get the AIC value and add it to the vector for all AIC values
	TheAIC=AIC(TheModel)
	AICs=c(AICs,TheAIC)
	
	# get the AIC value and add it to the vector for all AIC values
	Deviances=c(Deviances,TheModel$deviance)
	NullDeviances=c(NullDeviances,TheModel$null.deviance)
	
	# record the degrees of freedom (number of parameters) for each independent variable
	EDFs=c(EDFs,sum(TheModel$edf))
	
	Count=Count+1 # move to the next gamma value
}
PercentDeviance=Deviances/NullDeviances # find the percent deviance from the data

# plot the statistics from our MC run
par(mfrow=c(2,2)) # allow 2x2 graphs
plot(Gammas)
plot(AICs)
plot(PercentDeviance)
plot(EDFs) # plot the total effective degrees of freedom (number of parameters)