Last active
January 11, 2018 13:44
-
-
Save fawda123/4476610 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#import relevant packages | |
library(gstat) | |
library(sp) | |
library(nlme) | |
#create simulated data, lat and lon from uniformly distributed variable, exp1 and exp2 from random normal | |
set.seed(2) | |
samp.sz<-400 | |
lat<-runif(samp.sz,-4,4) | |
lon<-runif(samp.sz,-4,4) | |
exp1<-rnorm(samp.sz) | |
exp2<-rnorm(samp.sz) | |
#resp is linear combination of explanatory variables plus an error term that is normally distributed | |
resp<-1+4*exp1-3*exp2-2*lat+rnorm(samp.sz) | |
#pair plots of variables | |
plot(data.frame(resp,exp1,exp2,lat)) | |
#correlation between variables | |
cor(data.frame(resp,exp1,exp2,lat)) | |
#get mods, check parameter values with actual, but we don't know spatial correlation with lat | |
mod<-lm(resp~exp1+exp2) | |
coefficients(summary(mod)) | |
#function for testing parameter values against actual | |
t_test<-function(x.bar,mu,se,deg.f,return.t=F){ | |
if(return.t) return(round((x.bar-mu)/se,3)) | |
pt(abs((x.bar-mu)/se),deg.f,lower.tail=F)*2 | |
} | |
#for first explanatory variable | |
t_test(3.8314, 4, 0.2534, 397) | |
#evidence of spatial correlation using bubble plot | |
dat<-data.frame(lon,lat,resids=resid(mod)) | |
coordinates(dat)<-c('lon','lat') | |
bubble(dat,zcol='resids') | |
#evidence of spatial correlation using variogram | |
var.mod<-variogram(resids~1,data=dat,alpha=c(0,45,90,135)) | |
plot(var.mod) | |
#refit model with correlation structure | |
mod.cor<-gls(resp~exp1+exp2,correlation=corGaus(form=~lat,nugget=TRUE)) | |
summary(mod.cor) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment