# R-skript for Bertoldi et al. (2010). Topographical and ecohydrological controls on land surface temperature in an alpine catchment # done by Andrea Bou-Vinals, Johannes Oehm, Daniela Straube rm(list=ls()) #remove everything setwd("I:\\Eigene Dateien\\Adv_Statistics2010\\presentation\\") install.packages ("foreign") library(foreign) bertoldi <- read.spss(file= "I:\\Eigene Dateien\\Adv_Statistics2010\\presentation\\datenpool_land_surface_temperature_lst_bertoldi.sav",to.data.frame = TRUE) bertoldi.withoutna<-na.omit(bertoldi) attach(bertoldi.withoutna) bertoldi.withoutna ### OLS 1 für LST g #### ols1 <- lm(LST_G_new ~ sin_inc_new+LU_6_NEW+DTM_new+soil_moisture_new+LU_1_NEW+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(ols1) install.packages ("QuantPsyc") library(QuantPsyc) beta1<-lm.beta(ols1) beta1 par(mfrow = c(2,2)) plot(ols1) #Robust regression #da dieses package einen falschen link aufweise, hab ich es auf der homepage (http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/2.13/)downgeloaded und über Pakete/Lade Paket direkt aufgerufen. library(MASS) rlm1 <- rlm(LST_G_new ~ sin_inc_new+LU_6_NEW+DTM_new+soil_moisture_new+LU_1_NEW+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(rlm1) # Test of multicollinearity install.packages ("car") library(car) vif(ols1) tolerance <- 1/vif(ols1) tolerance #Part correlation # variable sinus of sun incidence angle rsquare_fullmodel<-(1-var(ols1$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) ols1_wo1<-lm(LST_G_new ~ LU_6_NEW+DTM_new+soil_moisture_new+LU_1_NEW+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(ols1_wo1) ols1_o1<-lm(LST_G_new ~ sin_inc_new) summary(ols1_o1) rsquare_ab<-(1-var(ols1_o1$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo1$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # variable elevation ols1_wo2<-lm(LST_G_new ~ LU_6_NEW+sin_inc_new+soil_moisture_new+LU_1_NEW+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(ols1_wo2) ols1_o2<-lm(LST_G_new ~ DTM_new) summary(ols1_o2) rsquare_ab<-(1-var(ols1_o2$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo2$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable Glaciers ols1_wo3<-lm(LST_G_new ~soil_moisture_new +sin_inc_new+DTM_new+LU_1_NEW+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(ols1_wo3) ols1_o3<-lm(LST_G_new ~ LU_6_NEW) summary(ols1_o3) rsquare_ab<-(1-var(ols1_o3$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo3$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable Soil moisture ols1_wo4<-lm(LST_G_new ~ LU_6_NEW+sin_inc_new+DTM_new+LU_1_NEW+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(ols1_wo4) ols1_o4<-lm(LST_G_new ~ soil_moisture_new) summary(ols1_o4) rsquare_ab<-(1-var(ols1_o4$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo4$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable forest ols1_wo5<-lm(LST_G_new ~ LU_6_NEW+sin_inc_new+DTM_new+soil_moisture_new+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(ols1_wo5) ols1_o5<-lm(LST_G_new ~ LU_1_NEW) summary(ols1_o5) rsquare_ab<-(1-var(ols1_o5$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo5$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable managed grassland ols1_wo6<-lm(LST_G_new ~ LU_6_NEW+sin_inc_new+DTM_new+soil_moisture_new+LU_1_NEW+LU_3_NEW+LU_4_NEW) summary(ols1_wo6) ols1_o6<-lm(LST_G_new ~ LU_2_NEW) summary(ols1_o6) rsquare_ab<-(1-var(ols1_o6$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo6$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable natural/abandoned grassland ols1_wo7<-lm(LST_G_new ~ LU_6_NEW+sin_inc_new+DTM_new+soil_moisture_new+LU_1_NEW+LU_2_NEW+LU_4_NEW) summary(ols1_wo7) ols1_o7<-lm(LST_G_new ~ LU_3_NEW) summary(ols1_o7) rsquare_ab<-(1-var(ols1_o7$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo7$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable settlement ols1_wo8<-lm(LST_G_new ~ LU_6_NEW+sin_inc_new+DTM_new+soil_moisture_new+LU_1_NEW+LU_2_NEW+LU_3_NEW) summary(ols1_wo8) ols1_o8<-lm(LST_G_new ~LU_4_NEW ) summary(ols1_o8) rsquare_ab<-(1-var(ols1_o8$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo8$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_G_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new #### - ols 2 für LSTL #### ols2 <- lm(LST_L_new ~ sin_inc_new+LU_6_NEW+DTM_new+LU_3_NEW+LU_1_NEW+LU_2_NEW+LU_4_NEW+soil_moisture_new) summary(ols2) beta2<-lm.beta(ols2) beta2 par(mfrow = c(2,2)) plot(ols2) #Robust regression rlm2 <- rlm(LST_L_new ~ sin_inc_new+LU_6_NEW+DTM_new+LU_3_NEW+LU_1_NEW+LU_2_NEW+LU_4_NEW+soil_moisture_new) summary(rlm2) beta2r<-lm.beta(rlm2) beta2r vif(ols2) tolerance <- 1/vif(ols2) tolerance #Part correlation # variable sinus of sun incidence angle rsquare_fullmodel2<-(1-var(ols2$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) ols2_wo1<-lm(LST_L_new ~ LU_6_NEW+DTM_new+soil_moisture_new+LU_1_NEW+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(ols2_wo1) ols2_o1<-lm(LST_L_new ~ sin_inc_new) summary(ols2_o1) rsquare_ab<-(1-var(ols2_o1$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols2_wo1$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel2) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # variable elevation #### ols1_wo2<-lm(LST_L_new ~ LU_6_NEW+sin_inc_new+soil_moisture_new+LU_1_NEW+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(ols1_wo2) ols1_o2<-lm(LST_L_new ~ DTM_new) summary(ols1_o2) rsquare_ab<-(1-var(ols1_o2$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo2$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel2) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable Glaciers ols1_wo3<-lm(LST_L_new ~soil_moisture_new +sin_inc_new+DTM_new+LU_1_NEW+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(ols1_wo3) ols1_o3<-lm(LST_L_new ~ LU_6_NEW) summary(ols1_o3) rsquare_ab<-(1-var(ols1_o3$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo3$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel2) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable Soil moisture ols1_wo4<-lm(LST_L_new ~ LU_6_NEW+sin_inc_new+DTM_new+LU_1_NEW+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(ols1_wo4) ols1_o4<-lm(LST_L_new ~ soil_moisture_new) summary(ols1_o4) rsquare_ab<-(1-var(ols1_o4$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo4$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel2) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable forest ols1_wo5<-lm(LST_L_new ~ LU_6_NEW+sin_inc_new+DTM_new+soil_moisture_new+LU_2_NEW+LU_3_NEW+LU_4_NEW) summary(ols1_wo5) ols1_o5<-lm(LST_L_new ~ LU_1_NEW) summary(ols1_o5) rsquare_ab<-(1-var(ols1_o5$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo5$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel2) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable managed grassland ols1_wo6<-lm(LST_L_new ~ LU_6_NEW+sin_inc_new+DTM_new+soil_moisture_new+LU_1_NEW+LU_3_NEW+LU_4_NEW) summary(ols1_wo6) ols1_o6<-lm(LST_L_new ~ LU_2_NEW) summary(ols1_o6) rsquare_ab<-(1-var(ols1_o6$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo6$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel2) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable natural/abandoned grassland ols1_wo7<-lm(LST_L_new ~ LU_6_NEW+sin_inc_new+DTM_new+soil_moisture_new+LU_1_NEW+LU_2_NEW+LU_4_NEW) summary(ols1_wo7) ols1_o7<-lm(LST_L_new ~ LU_3_NEW) summary(ols1_o7) rsquare_ab<-(1-var(ols1_o7$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo7$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel2) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new # Variable settlement ols1_wo8<-lm(LST_L_new ~ LU_6_NEW+sin_inc_new+DTM_new+soil_moisture_new+LU_1_NEW+LU_2_NEW+LU_3_NEW) summary(ols1_wo8) ols1_o8<-lm(LST_L_new ~LU_4_NEW ) summary(ols1_o8) rsquare_ab<-(1-var(ols1_o8$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) rsquare_bc<-(1-var(ols1_wo8$resid,na.rm = TRUE)/var(bertoldi.withoutna$LST_L_new,na.rm = TRUE)) partcorr_sin_inc_new_squared<-rsquare_ab-(rsquare_ab+rsquare_bc-rsquare_fullmodel2) partcorr_sin_inc_new<-partcorr_sin_inc_new_squared^0.5 partcorr_sin_inc_new ### subset #### rm(list=ls()) subset1<- subset(bertoldi.withoutna, LST_L_new>10) attach(subset1) # no values for variable LU_6_new (glacier), therefore this variable ist not included in lm-analysis ols_sub1 <- lm(LST_L_new ~ sin_inc_new+DTM_new+LU_3_NEW+LU_1_NEW+LU_2_NEW+LU_4_NEW+soil_moisture_new) summary(ols_sub1) par(mfrow = c(2,2)) plot(ols_sub1) #Robust regression rlm_sub1 <- rlm(LST_L_new ~ sin_inc_new+DTM_new+LU_3_NEW+LU_1_NEW+LU_2_NEW+LU_4_NEW+soil_moisture_new) summary(rlm_sub1) #### subset2<- subset(bertoldi.withoutna, LST_L_new<10) attach(subset2) ols_sub2 <- lm(LST_L_new ~ sin_inc_new+LU_6_NEW+DTM_new+LU_3_NEW+LU_1_NEW+LU_2_NEW+LU_4_NEW+soil_moisture_new) summary(ols_sub2) par(mfrow = c(2,2)) plot(ols_sub2) #Robust regression rlm_sub2 <- rlm(LST_L_new ~ sin_inc_new+LU_6_NEW+DTM_new+LU_3_NEW+LU_1_NEW+LU_2_NEW+soil_moisture_new) summary(rlm_sub2)