Koden finns även tillgänglig på GitHub via https://tinyurl.com/prediktionR.
OBS: Notera att koden inte är exekverbar utan tillgång till de tretton GIS-projekt arbetet är baserat på, vilka överstiger 10 GB och således inte kommer bifogas. I och med att koden ursprungligen endast var avsett för detta arbete är koden dåligt formaterad och sporadiskt kommenterad. Merparten av koden bör dock kunna appliceras i andra prediktionsmodeller så länge sökvägar till filerna korrigeras.
install.packages("devtools")
install.packages("lme4", repos="http://cran.rstudio.com/",type = "binary", dependencies=TRUE) install.packages("nlme", repos="http://cran.rstudio.com/",type = "binary", dependencies=TRUE) packageurl <- "https://cran.r-project.org/src/contrib/Archive/pbkrtest/pbkrtest_0.4-4.tar.gz"
preddata <- readOGR('anvand/final.shp') #undvik att kora upprepade ganger summary(preddata)
unique(preddata$JG2_TX) vegrelfreq <- vegfreq / nrow(testdata) bergrelfreq <- bergfreq / nrow(testdata) jordrelfreq <- jordfreq / nrow(testdata) vegfreq
testdata$amfi[which (testdata$Bergart!="Amfibolit")] <- 0 testdata$amfi[which (testdata$Bergart=="Amfibolit")] <- 1
testdata$gabb_dior <- NA
testdata$gabb_dior[which (testdata$Bergart!="Gabbroid-dioritoid")] <- 0 testdata$gabb_dior[which (testdata$Bergart=="Gabbroid-dioritoid")] <- 1
testdata$gran <- NA
testdata$gran[which (testdata$Bergart!="Granit")] <- 0 testdata$gran[which (testdata$Bergart=="Granit")] <- 1
testdata$gran_gran <- NA
testdata$gran_gran[which (testdata$Bergart!="Granodiorit-granit")] <- 0 testdata$gran_gran[which (testdata$Bergart=="Granodiorit-granit")] <- 1
testdata$kalk <- NA
testdata$kalk[which (testdata$Bergart!="Kalksten")] <- 0 testdata$kalk[which (testdata$Bergart=="Kalksten")] <- 1
testdata$kvar <- NA
testdata$kvar[which (testdata$Bergart!="Kvartsarenit")] <- 0 testdata$kvar[which (testdata$Bergart=="Kvartsarenit")] <- 1
testdata$para <- NA
testdata$para[which (testdata$Bergart!="Paragnejs")] <- 0 testdata$para[which (testdata$Bergart=="Paragnejs")] <- 1
testdata$sands <- NA
testdata$sands[which (testdata$Bergart!="Sandsten")] <- 0 testdata$sands[which (testdata$Bergart=="Sandsten")] <- 1
testdata$skif <- NA
testdata$skif[which (testdata$Bergart!="Skiffer")] <- 0 testdata$skif[which (testdata$Bergart=="Skiffer")] <- 1
testdata$syen_gran <- NA
testdata$syen_gran[which (testdata$Bergart!="Syenitoid-granit")] <- 0 testdata$syen_gran[which (testdata$Bergart=="Syenitoid-granit")] <- 1
testdata$tona_gran <- NA
testdata$tona_gran[which (testdata$Bergart!="Tonalit-granodiorit")] <- 0 testdata$tona_gran[which (testdata$Bergart=="Tonalit-granodiorit")] <- 1
testdata$vacka <- NA
testdata$vacka[which (testdata$Bergart!="Vacka")] <- 0 testdata$vacka[which (testdata$Bergart=="Vacka")] <- 1
preddata$amfi <- NA
preddata$amfi[which (preddata$BERGART_TX!="Amfibolit")] <- 0 preddata$amfi[which (preddata$BERGART_TX=="Amfibolit")] <- 1
preddata$gabb_dior <- NA
preddata$gabb_dior[which (preddata$BERGART_TX!="Gabbroid-dioritoid")] <- 0 preddata$gabb_dior[which (preddata$BERGART_TX=="Gabbroid-dioritoid")] <- 1
preddata$gran <- NA
preddata$gran[which (preddata$BERGART_TX!="Granit")] <- 0 preddata$gran[which (preddata$BERGART_TX=="Granit")] <- 1
preddata$gran_gran <- NA
preddata$gran_gran[which (preddata$BERGART_TX!="Granodiorit-granit")] <- 0 preddata$gran_gran[which (preddata$BERGART_TX=="Granodiorit-granit")] <- 1
preddata$kalk <- NA
preddata$kalk[which (preddata$BERGART_TX!="Kalksten")] <- 0 preddata$kalk[which (preddata$BERGART_TX=="Kalksten")] <- 1
preddata$kvar <- NA
preddata$kvar[which (preddata$BERGART_TX!="Kvartsarenit")] <- 0 preddata$kvar[which (preddata$BERGART_TX=="Kvartsarenit")] <- 1
preddata$para <- NA
preddata$para[which (preddata$BERGART_TX!="Paragnejs")] <- 0 preddata$para[which (preddata$BERGART_TX=="Paragnejs")] <- 1
preddata$sands <- NA
preddata$sands[which (preddata$BERGART_TX!="Sandsten")] <- 0 preddata$sands[which (preddata$BERGART_TX=="Sandsten")] <- 1
preddata$skif <- NA
preddata$skif[which (preddata$BERGART_TX!="Skiffer")] <- 0 preddata$skif[which (preddata$BERGART_TX=="Skiffer")] <- 1
preddata$syen_gran <- NA
preddata$syen_gran[which (preddata$BERGART_TX!="Syenitoid-granit")] <- 0 preddata$syen_gran[which (preddata$BERGART_TX=="Syenitoid-granit")] <- 1
preddata$tona_gran <- NA
preddata$tona_gran[which (preddata$BERGART_TX!="Tonalit-granodiorit")] <- 0 preddata$tona_gran[which (preddata$BERGART_TX=="Tonalit-granodiorit")] <- 1
preddata$vacka <- NA
preddata$vacka[which (preddata$BERGART_TX!="Vacka")] <- 0 preddata$vacka[which (preddata$BERGART_TX=="Vacka")] <- 1
#binara jordarter
testdata$berg <- NA
testdata$berg[which (testdata$Jordart!="Berg")] <- 0 testdata$berg[which (testdata$Jordart=="Berg")] <- 1
testdata$isalv <- NA
testdata$isalv[which (testdata$Jordart!="Isälvssediment")] <- 0 testdata$isalv[which (testdata$Jordart=="Isälvssediment")] <- 1
testdata$lera_silt<- NA
testdata$lera_silt[which (testdata$Jordart!="Lera--silt")] <- 0 testdata$lera_silt[which (testdata$Jordart=="Lera--silt")] <- 1
testdata$moran <- NA
testdata$moran[which (testdata$Jordart!="Morän")] <- 0 testdata$moran[which (testdata$Jordart=="Morän")] <- 1
testdata$torv<- NA
testdata$torv[which (testdata$Jordart!="Torv")] <- 0 testdata$torv[which (testdata$Jordart=="Torv")] <- 1
testdata$vattenjord<- NA
testdata$vattenjord[which (testdata$Jordart!="Vatten")] <- 0
testdata$vattenjord[which (testdata$Jordart=="Vatten")] <- 1
testdata$alvs_sand<- NA
testdata$alvs_sand[which (testdata$Jordart!="Älvsediment, sand")] <- 0 testdata$alvs_sand[which (testdata$Jordart=="Älvsediment, sand")] <- 1
unique(preddata$JG2_TX) preddata$berg <- NA
preddata$berg[which (preddata$JG2_TX!="Berg")] <- 0 preddata$berg[which (preddata$JG2_TX=="Berg")] <- 1
preddata$isalv <- NA
preddata$isalv[which (preddata$JG2_TX!="Is?lvssediment")] <- 0 preddata$isalv[which (preddata$JG2_TX=="Is?lvssediment")] <- 1
preddata$lera_silt<- NA
preddata$lera_silt[which (preddata$JG2_TX!="Lera--silt")] <- 0 preddata$lera_silt[which (preddata$JG2_TX=="Lera--silt")] <- 1
preddata$moran <- NA
preddata$moran[which (preddata$JG2_TX!="Mor?n")] <- 0 preddata$moran[which (preddata$JG2_TX=="Mor?n")] <- 1
preddata$torv<- NA
preddata$torv[which (preddata$JG2_TX!="Torv")] <- 0 preddata$torv[which (preddata$JG2_TX=="Torv")] <- 1
preddata$vattenjord<- NA
preddata$vattenjord[which (preddata$JG2_TX!="Vatten")] <- 0 preddata$vattenjord[which (preddata$JG2_TX=="Vatten")] <- 1
preddata$alvs_sand<- NA
preddata$alvs_sand[which (preddata$JG2_TX!="Älvsediment, sand")] <- 0 preddata$alvs_sand[which (preddata$JG2_TX=="Älvsediment, sand")] <- 1
#binar vegetation vegfreq
testdata$barr_lav<- NA
testdata$barr_lav[which (testdata$Vegetationstyp!="Barrskog, lavristyp")] <- 0 testdata$barr_lav[which (testdata$Vegetationstyp=="Barrskog, lavristyp")] <- 1
testdata$myr_karr<- NA
testdata$myr_karr[which (testdata$Vegetationstyp!="Barrskogsmyr")] <- 0 testdata$myr_karr[which (testdata$Vegetationstyp=="Barrskogsmyr")] <- 1
testdata$expl<- NA
testdata$expl[which (testdata$Vegetationstyp!="Exploaterad mark")] <- 0 testdata$expl[which (testdata$Vegetationstyp=="Exploaterad mark")] <- 1
testdata$myr_karr[which (testdata$Vegetationstyp=="Fastmattemyr, halvgräsvariant")] <- 1
testdata$barr_fukt<- NA
testdata$barr_fukt[which (testdata$Vegetationstyp!="Fuktig barrskog")] <- 0 testdata$barr_fukt[which (testdata$Vegetationstyp=="Fuktig barrskog")] <- 1
testdata$myr_karr[which (testdata$Vegetationstyp=="Högstarr-sumpkärr")] <- 1
testdata$kult<- NA
testdata$kult[which (testdata$Vegetationstyp!="Kulturmark")] <- 0 testdata$kult[which (testdata$Vegetationstyp=="Kulturmark")] <- 1
testdata$myr_karr[which (testdata$Vegetationstyp=="Lösbottenkärr, starr-örtvariant")] <- 1
testdata$myr_karr[which (testdata$Vegetationstyp=="Lösbottenmyr, halvgräsvariant")] <- 1
testdata$myr_karr[which (testdata$Vegetationstyp==" Mjukmattemyr, halvgräs-vitmossvariant")] <- 1
testdata$barr_moss<- NA
testdata$barr_moss[which (testdata$Vegetationstyp!="Mossmarksbarrskog")] <- 0 testdata$barr_moss[which (testdata$Vegetationstyp=="Mossmarksbarrskog")] <- 1
testdata$myr_karr[which (testdata$Vegetationstyp=="Ristuvemyr")] <- 1
testdata$strand<- NA
testdata$strand[which (testdata$Vegetationstyp!="Sötvattensstrandäng (nordlig), sedimentationsbetingad")] <- 0 testdata$strand[which (testdata$Vegetationstyp=="Sötvattensstrandäng (nordlig), sedimentationsbetingad")] <- 1
testdata$barr_torr<- NA
testdata$barr_torr[which (testdata$Vegetationstyp!="Torr-frisk barrskog")] <- 0 testdata$barr_torr[which (testdata$Vegetationstyp=="Torr-frisk barrskog")] <- 1
testdata$barr_vat<- NA
testdata$barr_vat[which (testdata$Vegetationstyp!="Våt barrskog")] <- 0 testdata$barr_vat[which (testdata$Vegetationstyp=="Våt barrskog")] <- 1
testdata$vattenveg<- NA
testdata$vattenveg[which (testdata$Vegetationstyp!="Öppet vatten")] <- 0 testdata$vattenveg[which (testdata$Vegetationstyp=="Öppet vatten")] <- 1
unique(preddata$VEGETATION)
preddata$barr_lav<- NA
preddata$barr_lav[which (preddata$VEGETATION!="Barrskog, lavristyp")] <- 0 preddata$barr_lav[which (preddata$VEGETATION=="Barrskog, lavristyp")] <- 1
preddata$myr_karr<- NA
preddata$myr_karr[which (preddata$VEGETATION!="Barrskogsmyr")] <- 0 preddata$myr_karr[which (preddata$VEGETATION=="Barrskogsmyr")] <- 1
preddata$expl<- NA
preddata$expl[which (preddata$VEGETATION!="Exploaterad mark")] <- 0 preddata$expl[which (preddata$VEGETATION=="Exploaterad mark")] <- 1
preddata$myr_karr[which (preddata$VEGETATION=="Fastmattemyr, halvgr?svariant")] <- 1
preddata$barr_fukt<- NA
preddata$barr_fukt[which (preddata$VEGETATION!="Fuktig barrskog")] <- 0 preddata$barr_fukt[which (preddata$VEGETATION=="Fuktig barrskog")] <- 1
preddata$myr_karr[which (preddata$VEGETATION=="Högstarr-sumpkärr")] <- 1
preddata$kult<- NA
preddata$kult[which (preddata$VEGETATION!="Kulturmark")] <- 0 preddata$kult[which (preddata$VEGETATION=="Kulturmark")] <- 1
preddata$myr_karr[which (preddata$VEGETATION=="L?sbottenk?rr, starr-?rtvariant")] <- 1
preddata$myr_karr[which (preddata$VEGETATION=="Lösbottenmyr, halvgräsvariant")] <- 1
preddata$myr_karr[which (preddata$VEGETATION==" Mjukmattemyr, halvgr?s-vitmossvariant")] <- 1
preddata$barr_moss<- NA
preddata$barr_moss[which (preddata$VEGETATION!="Mossmarksbarrskog")] <- 0 preddata$barr_moss[which (preddata$VEGETATION=="Mossmarksbarrskog")] <- 1
preddata$myr_karr[which (preddata$VEGETATION=="Ristuvemyr")] <- 1
preddata$strand<- NA
preddata$strand[which (preddata$VEGETATION!="Sötvattensstrandäng (nordlig), sedimentationsbetingad")]
<- 0
preddata$strand[which (preddata$VEGETATION=="Sötvattensstrandäng (nordlig), sedimentationsbetingad")]
<- 1
preddata$barr_torr<- NA
preddata$barr_torr[which (preddata$VEGETATION!="Torr-frisk barrskog")] <- 0 preddata$barr_torr[which (preddata$VEGETATION=="Torr-frisk barrskog")] <- 1
preddata$barr_vat<- NA
preddata$barr_vat[which (preddata$VEGETATION!="V?t barrskog")] <- 0 preddata$barr_vat[which (preddata$VEGETATION=="V?t barrskog")] <- 1
preddata$vattenveg<- NA
preddata$vattenveg[which (preddata$VEGETATION!="?ppet vatten")] <- 0 preddata$vattenveg[which (preddata$VEGETATION=="?ppet vatten")] <- 1
#slut pa binar kodning
sum(is.na(testdata$Jordart)) sum(is.na(testdata$Bergart)) sum(is.na(testdata$Vegetationstyp))
#spplot(preddata, "VEGETATION")
lmModAvstand <- lm(testdata$Bobin ~ testdata$Avstand, data=trainingData) # build the model lmModSlope <- lm(testdata$Bobin ~ testdata$Slopeint, data=trainingData)
lmModView <- lm(testdata$Bobin ~ testdata$View, data=trainingData) lmModHojd <- lm(testdata$Bobin ~ testdata$Omkodad, data=trainingData) lmModamfi <- lm(testdata$Bobin ~ testdata$amfi, data=trainingData)
lmModgabb_dior <- lm(testdata$Bobin ~ testdata$gabb_dior, data=trainingData) lmModgran <- lm(testdata$Bobin ~ testdata$gran, data=trainingData)
lmModgran_gran <- lm(testdata$Bobin ~ testdata$gran_gran, data=trainingData) lmModkalk <- lm(testdata$Bobin ~ testdata$kalk, data=trainingData)
lmModkvar <- lm(testdata$Bobin ~ testdata$kvar, data=trainingData) lmModpara <- lm(testdata$Bobin ~ testdata$para, data=trainingData) lmModsands <- lm(testdata$Bobin ~ testdata$sands, data=trainingData) lmModskif <- lm(testdata$Bobin ~ testdata$skif, data=trainingData)
lmModsyen_gran <- lm(testdata$Bobin ~ testdata$syen_gran, data=trainingData) lmModtona_gran <- lm(testdata$Bobin ~ testdata$tona_gran, data=trainingData) lmModvacka <- lm(testdata$Bobin ~ testdata$vacka, data=trainingData) lmModberg <- lm(testdata$Bobin ~ testdata$berg, data=trainingData) lmModisalv <- lm(testdata$Bobin ~ testdata$isalv, data=trainingData) lmModlera_silt <- lm(testdata$Bobin ~ testdata$lera_silt, data=trainingData) lmModmoran <- lm(testdata$Bobin ~ testdata$moran, data=trainingData) lmModtorv <- lm(testdata$Bobin ~ testdata$torv, data=trainingData)
lmModvattenjord <- lm(testdata$Bobin ~ testdata$vattenjord, data=trainingData) lmModalvs_sand <- lm(testdata$Bobin ~ testdata$alvs_sand, data=trainingData) lmModbarr_lav <- lm(testdata$Bobin ~ testdata$barr_lav, data=trainingData) lmModexpl <- lm(testdata$Bobin ~ testdata$expl, data=trainingData)
lmModmyr_karr <- lm(testdata$Bobin ~ testdata$myr_karr, data=trainingData) lmModbarr_fukt <- lm(testdata$Bobin ~ testdata$barr_fukt, data=trainingData) lmModkult <- lm(testdata$Bobin ~ testdata$kult, data=trainingData)
lmModbarr_moss <- lm(testdata$Bobin ~ testdata$barr_moss, data=trainingData) lmModstrand <- lm(testdata$Bobin ~ testdata$strand, data=trainingData)
lmModbarr_torr <- lm(testdata$Bobin ~ testdata$barr_torr, data=trainingData) lmModbarr_vat <- lm(testdata$Bobin ~ testdata$barr_vat, data=trainingData) lmModvattenjord <- lm(testdata$Bobin ~ testdata$vattenjord, data=trainingData)
preddata$distPred <- predict(fit, preddata) # predict distance
fit <- lm(Bobin~amfi + gabb_dior + gran+gran_gran+kalk+kvar+para +sands +skif +syen_gran +tona_gran +vacka +berg +isalv +lera_silt +moran +torv +vattenjord +alvs_sand +barr_lav +expl +myr_karr +barr_fukt +kult +barr_moss +strand +barr_torr +barr_vat+vattenjord+Omkodad+View +Slopeint + Avstand , data=testdata)
fitexp <- NA
fitexp <- glm(Bobin~amfi + gabb_dior + gran+gran_gran+kalk+kvar+para +sands +skif +syen_gran +tona_gran +vacka +berg +isalv +lera_silt +moran +torv +vattenjord +alvs_sand +barr_lav +expl +myr_karr +barr_fukt +kult +barr_moss +strand +barr_torr +barr_vat+vattenjord* Omkodad* View * Slopeint * Avstand , data=testdata)
predfit <- NA
predfit$barr_fukt <- preddata$barr_fukt
writeOGR(preddata, dsn = "anvand/anvand.shp", layer = "test", driver = "ESRI Shapefile") summary(pred)
spplot(preddata, "pred")
spplot(preddata)
set.seed()
training.samples <- testdata$Bobin %>%
createDataPartition(p = 0.8, list = FALSE) train.data <- testdata[training.samples, ] test.data <- testdata[-training.samples, ]
model <- lm(Bobin~amfi + gabb_dior + gran+gran_gran+kalk+kvar+para +sands +skif +syen_gran +tona_gran +vacka +berg +isalv +lera_silt +moran +torv +vattenjord +alvs_sand +barr_lav +expl +myr_karr +barr_fukt +kult +barr_moss +strand +barr_torr +barr_vat+vattenjord+Omkodad+View +Slopeint + Avstand , data=testdata)
predictions <- model %>% predict(test.data) data.frame( R2 = R2(predictions, test.data$Bobin), RMSE = RMSE(predictions, test.data$Bobin), MAE = MAE(predictions, test.data$Bobin)) model <- NA
train.control <- trainControl(method = "LOOCV") model <- train(Bobin ~., data = testdata, method = "lm", trControl = train.control, na.action = NULL)
print(model)
train.control <- trainControl(method = "LOOCV") model <- train(Bobin ~., data = testdata, method = "lm", trControl = train.control, na.action = NULL)
print(model)
fit
outlierTest(fit)
qqPlot(fit, main="QQ Plot") leveragePlots(fit)
av.Plots(fit)
cutoff <- 4/((nrow(mtcars)-length(fit$coefficients)-2)) plot(fit, which=4, cook.levels=cutoff)
influencePlot(fit, id.method="identify", main="Influence Plot", sub="Circle size is proportial to Cook's Distance" )