This lab will introduce you to modeling presence/absence data with GLMs. This is also the first full-modeling lab. A key element of this lab is examining the response of the data vs. the predicted output of the model and how both relate to the predictor variables.
Note: Remember to include the glm2 library in the code below.
The function below will create synthetic presence/absence data for evaluating GLMs. Copy it into R now and compile them.
############################################################################
# Creates a data frame with Xs ranging from 1 to the number of entries
# and "Measures" with 1/2 set to 0 and 1/2 set to 1. The lower have of
# the X values have measures at 0.
# ProportionRandom - amount of uniform randomness to add
############################################################################
Categories1D_Random=function(NumEntries=10,ProportionRandom=0.4)
{
Range=NumEntries*ProportionRandom/2
Ys=as.vector(array(1:NumEntries))
Xs=as.vector(array(1:NumEntries))
Measures=as.vector(array(1:NumEntries))
for (Index in 1:NumEntries)
{
Xs[Index]=Index #runif(1,0,100)
Ys[Index]=Index #runif(1,0,100)
Threshold=0.5
Random=0
if (ProportionRandom!=0) Random=runif(1,-Range,Range)
if (Xs[Index]>NumEntries/2+Random) Measures[Index]=1
else Measures[Index]=0
}
TheDataFrame = data.frame(Ys, Xs, Measures)
}
The code below will create a synthetic data set for a logistic model using the function above. Try it now in R. Note that the data contains some randomness to make the values overlap a bit.
TheData=Categories1D_Random(100) # create a set of binary data plot(TheData$Xs,TheData$Measures) # plot the data
The code below will fit the model and plot the results.
TheModel=glm2(Measures~Xs,data=TheData,family=binomial) # logit (binary) summary(TheModel)
The output is similar to a linear model with the exception of a few new items.
Call:
glm2(formula = Measures ~ Xs, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.19011 -0.18968 0.02669 0.18966 2.25542
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.86152 1.78795 -4.397 1.10e-05 ***
Xs 0.15882 0.03521 4.511 6.45e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 138.589 on 99 degrees of freedom
Residual deviance: 41.339 on 98 degrees of freedom
AIC: 45.339
Number of Fisher Scoring iterations: 7
Note that at the bottom is the AIC. This and the model itself, are the two key pieces of information for model selection. The "Null deviance" is effectively the "-2ln(L)" portion of the AIC for a "Null" model and the "Residual deviance" is the same calculation for the fitted model (i.e. the residual deviance should be less than the Null deviance. The "Fisher Scoring" is the method used to find the maximum likelihood and is not critical for model selection.
I'm not sure but I think: The standard deviation of the ceofficient represents how well the coeficient has been estimated (http://tinyurl.com/j5ccvgo). The z value is the regression coeficient divided by the standard error. If the z-value is large (>2), then the coeficient matters (http://tinyurl.com/hxexbtg). The Pr(>|z|) is like a p-value.
Note that the plots provided are described below, however, these plots are intended for linear models and can be misleading for GLMs.
plot(TheModel)
We can obtain the residuals and histogram them as below:
residuals(TheModel)
hist(residuals(TheModel))
We can view the fitted values and plot them with:
fitted(TheModel)
plot(fitted(TheModel))
Note that there are two functions that will create a prediction for a GLM, "predict()" and "fitted()". Sources on the web indicate that predict() does not include the "link" function so it should not be compared with the original data. However, in my experience, predict() does overlay the original data properly so it must be applying the link function and fitted() will only predict against the original model values. If you want to predict against a new set of predictors (to make a map, for instance), you'll need to use predict. Note that the "type="response"" parameter is also required.
Xs=as.vector(array(1:100)) # create a vector with the full range of x values NewData=data.frame(Xs) # put this into a dataframe for prediction (names must match the model) ThePrediction=predict(TheModel,newdata=NewData,type="response") # create the prediction plot(TheData$Xs,TheData$Measures) # plot the original points points(Xs,ThePrediction,col=2,pch="+") # plot the prediction
Note that GLM is much pickier about it's inputs than most R functions. You may receive errors for having too few points or too many. Also, if it cannot find a reasonable result, it can also provide errors.
Start by plotting your prediction against the original data and you should have something like the following:
plot(TheData$Xs,TheData$Measures) # plot the data points(Xs,ThePrediction,col=2,pch="+") # plot the prediction
Note that the residuals increase as we move from the response values of 0 and 1 toward 0.5. This is inherent in a model where we are using a logistic equation to model a binary phenomenon. Now, run the model plots again to produce four graphs that can be used to interpret the model:
plot(TheModel)
The "residuals vs. fitted" plot shows the residuals on the y-axis and the fitted values (y values from the samples) on the x-axis. Another way to think of this is that the x-axis shows the predicted value along the trend line and the y-axis shows the difference between this value and the sample value (i.e. the residual).
The next graph shows a Normal Q-Q plot. This plot breaks the original data into "quantiles" or buckets where each bucket as about the same number of data points. Then the mean value of each bucket is plotted for the original data against what we would expect from the model. The dashed line shows the expected value.
The problem with the graph above is that it is for the Standardized Deviance of the Residuals. I find this graph to be incomprehensible. Please do not include graphs unless you can explain what they mean.
The function used above to create the synthetic data has two parameters. The first indicates the number of data points to generate and the second the amount of randomness or "noise" to put into the model. Try some different values for these parameters including 0 for the randomness. Note that the GLM function cannot always create a valid model.
Below are some of the other link functions provided by the glm2() function. Try them with your data and evaluate the AIC values against the response curves provided.
TheModel=glm2(Measures~Xs,data=TheData,family=gaussian) # identity (linear) TheModel=glm2(Measures~Xs,data=TheData,family=poisson) # log (counts) TheModel=glm2(Measures~Xs,data=TheData,family=binomial) # logit (binary) TheModel=glm2(Measures~Xs,data=TheData,family=Gamma) # inverse (distances)
Note that the AIC values are not comparable across the methods for GLM. They are only valid for comparison within one method and with the same data set. We'll chat more on this later.
The files below represent examples of creating models for the presence of Doug-Fir trees in California. The Doug-Fir data set can be converted to "Points with Binary Values" as explained in the Converting Data Types web page and then adding predictor layers values from rasters for the study area. The code below is an example to create a predicted output for Doug-Fir based on the new presence-absence value. Note that these are examples and are not going to work exactly as you may desire. Also, you'll need to modify these examples for your own covariates or your own response data.
#########################################################
# Read and clean up the data
library(glm2)
# read the data into an R dataframe
TheData <- read.csv("C:/Temp/DougFir_CA2.csv")
# Remove any NA (null) values, the data should already be loaded
TheData=na.omit(TheData)
# plot the covariate against the response
plot(TheData$AnnualPrecip,TheData$PresAbs)
# create the model
TheModel=glm2(PresAbs~AnnualPrecip,data=TheData,family=binomial) # logit (binary)
summary(TheModel)
# create a new data frame for the prediction
AnnualPrecip=as.vector(array(1:2000)) # create a vector with the full range of x values
NewData=data.frame(AnnualPrecip) # put this into a dataframe for prediction (names must match the model)
# create the prediction fector
ThePrediction=predict(TheModel,newdata=NewData,type="response") # create the prediction
# plot the original points and the predicted points
plot(TheData$AnnualPrecip,TheData$PresAbs) # plot the original points
points(AnnualPrecip,ThePrediction,col=2,pch="+") # plot the prediction
# predict based on the original model values (i.e. the entire data set)
ThePrediction=predict(TheModel,type="response") # create the prediction
# add the prediction to the existing dataframe
FinalData=data.frame(TheData,ThePrediction)
plot(FinalData$AnnualPrecip,FinalData$ThePrediction)
# Save the prediction to a new file
write.csv(FinalData, "C:/Temp/ModelOutput.csv")
cv.glm, Cross-validation for Generalized Linear Models (uses k-fold)
There are some excellent labs available for R from MSU. Go through Lab 4 up to "Examples with the Use of R". The data for this lab is available in Lab 1. It is in an "R" file but needs to be loaded as a table as shown in the lab and there are some data quality issues that need to be addressed (I've made notes in a custom copy of the lab). It is also worthwhile to go through Lab 1 as practice.
© Copyright 2018 HSU - All rights reserved.