#### Problem 3 - Homework 1 for CVEN6833 - Fitting a General linear model for January
#### Author: Alvaro Ossandon
#### CLEAR MEMORY
rm(list=ls())
#### Clear the R console
cat("\f")
#### Prevent Warnings
options(warn=-1)
#data required:
Month=1 # 1: Monthly precipitation of January; 2: Monthly precipitation of July
save=FALSE # True for saving figure as pdf file;
#### SOURCE LIBRARIES (suppress package loading messages) AND SET WORKING DIRECTORY
mainDir="C:/Users/DELL/Dropbox/Boulder/Courses/Advance data analysis/HW/HW1"
setwd(mainDir)
suppressPackageStartupMessages(source("Lib_HW1.R"))
source("Lib_HW1.R")
#### create a folder to save the figures
if (Month==1){
diraux=paste(mainDir,"/Montly_Precep_January",sep = "")
}else{
diraux=paste(mainDir,"/Montly_Precep_July",sep = "")
}
setwd(diraux)
subDir="GLM_FIT"
if (file.exists(subDir)){
setwd(file.path(diraux, subDir))
} else {
dir.create(file.path(diraux, subDir))
setwd(file.path(diraux, subDir))
}
## Read data as table (data frame)
test = read.table(paste("http://cires1.colorado.edu/~aslater/CVEN_6833/colo_monthly_precip_0",Month,".dat",sep=""))
#assign names for columns/variables:
names(test)=c("Lat","Long","Elev","Pm")
test1=test
#filtering the -999 values
non_values<- which(test[,4]<0)
test[non_values,4]=NaN
test=na.omit(test)
Pm = test[,4] #Dependent Variable
X = test[,1:3]# Independent variable.
lon = X[,1]
lat = X[,2]
elev = X[,3]
precip =Pm
# Basic Scatterplot Matrix
if (save==TRUE){
pdf("Simple Scatterplot Matrix.pdf") # save figure
pairs(~lon+lat+elev+precip,data=test,
main="Simple Scatterplot Matrix")
dev.off()
} else{
pairs(~lon+lat+elev+precip,data=test,
main="Simple Scatterplot Matrix")
}
################################################
############## ii. Fit a GLM ###############
################################################
##Gamma and Gaussian distribution are tested
bestmod_gamma=GLM_fit(Pm,X,"gamma",Month)
## [1] "Results of PRESS for bestfit GLM"
## glm_press
## log 155.46225
## inverse 200000.00000
## identity 69.94976
## [1] 7
bestmod_gau=GLM_fit(Pm,X,"gaussian",Month)
## [1] "Results of PRESS for bestfit GLM"
## glm_press
## identity 265941.2
## [1] 7
# Test using PRESS objective function
if (bestmod_gamma$Call$PRESS<bestmod_gau$Call$PRESS){
bestmod=bestmod_gamma
bestmod$Call$Press=PRESS(bestmod)
} else{
bestmod=bestmod_gau
}
summary(bestmod)
##
## Call:
## glm(formula = Pm ~ ., family = Gamma(link = links[which.min(glm_press)]),
## data = xx)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.74780 -0.34667 -0.06703 0.15992 1.84470
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.213922 47.512389 0.404 0.686
## Lat 0.855067 0.616999 1.386 0.166
## Long 0.562926 0.398258 1.413 0.158
## Elev 0.021316 0.001446 14.741 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1890912)
##
## Null deviance: 140.963 on 448 degrees of freedom
## Residual deviance: 69.092 on 445 degrees of freedom
## AIC: 3660.8
##
## Number of Fisher Scoring iterations: 11
bestmod$Call$Press
## [1] 69.94976
Pmhat=bestmod$fitted.values
if (save==TRUE){
pdf("Precipitation_Scatterplot.pdf", width = 8, height = 6) # save figure
# Observed versus estimates
par(mfrow=c(1,1))
lim=range(Pm,Pmhat)
plot(Pm,Pmhat,xlab="Actual Precipitation",ylab="Estimated Precipitation",main="Actual vs Estimated Precipitation", xlim=lim, ylim=lim)
abline(a=0,b=1)
dev.off()
}else{
# Observed versus estimates
par(mfrow=c(1,1))
lim=range(Pm,Pmhat)
plot(Pm,Pmhat,xlab="Actual Precipitation",ylab="Estimated Precipitation",main="Actual vs Estimated Precipitation", xlim=lim, ylim=lim)
abline(a=0,b=1)
}
###################################################################
######## III. Anova and model diagnostics of your best model
###################################################################
### ANOVA:
anova(bestmod)
## Analysis of Deviance Table
##
## Model: Gamma, link: links[which.min(glm_press)]
##
## Response: Pm
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 448 140.963
## Lat 1 0.144 447 140.819
## Long 1 34.355 446 106.465
## Elev 1 37.373 445 69.092
### MODEL DIAGNOSTICS:
#yy = precip_dt[,4] # true value of y based on data
Pmest = Pmhat # model's predicted values of Pm
nvar = 2 # number of variables
mod_diagnostics(Pm, Pmest, nvar,save)
###################################################################
#### IV. Compute cross-validated and fitted estmiates at each ####
#### observation. Plot them against the observed values. ####
###################################################################
loocv_func(bestmod,X,Pm,save)
###################################################################
#### V. Compute cross-validated and fitted estmiates at each ####
#### observation. Plot them against the observed values. ####
###################################################################
Drop_10_pred(bestmod,X,Pm,save)
###################################################################
## VI. Spatially map the model estimates and the standard error ####
###################################################################
# read the topography map
x1 = read.table("http://cires1.colorado.edu/~aslater/CVEN_6833/colo_dem.dat")
names(x1)=c("Lat","Long","Elev")
ypred =predict.glm(bestmod, newdata=x1[,1:3], se.fit=T,type="response")
#Quilt_plotting(x1[,2],x1[,1],ypred)
plot_surface(x1[,2],x1[,1],ypred,X[,2],X[,1],Pm,save)
