##################################################################################
#
#		Traits prediction preparation - Sampling Sites - emprirical
#
#################################################################################

{
  library(Hmsc)
  library(tidyr)
  library(dplyr)
  library(abind)

}



## Trait info
#........................................................................................
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/2.ObsInfo")
Traits = read.csv('Traits_Moths.csv', header=T)
head(Traits)

Traits.sub= Traits %>%
  dplyr::select(Genus.and.species, female.size..average.wingspan..mm.,
                HostPlantForm)

Traits.sub$HostPlantForm = as.factor(Traits.sub$HostPlantForm)
levels(Traits.sub$HostPlantForm) <- c('herbaceous', 'graminoids', 'lichens-algae', 'Other', 'tree-stemmed', 
                                      'decideoustrees-shrub', 'Conifer-shrubs')

head(Traits.sub)
colnames(Traits.sub)[1] = 'species'
colnames(Traits.sub)[2] = 'WingspanF'
Traits.sub$WingGP = NA
Traits.sub$WingGP[Traits.sub$WingspanF <28] = 'W<28'
Traits.sub$WingGP[Traits.sub$WingspanF > 28 & Traits.sub$WingspanF < 34] = '28<W<34'
Traits.sub$WingGP[Traits.sub$WingspanF > 34] = 'W>34'


# cluster info
#........................................................................................
load('MothclusterFIN.RDATA')
dim(Moth_knn)  # 109    13
head(Moth_knn)
Moth_knn2 = Moth_knn[, 7:8]
colnames(Moth_knn2)[2] = 'sample'




## Results from predictions
#........................................................................................
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/4.Predictions/3.Predictions_Sites")

load('predABeveSIT_Y1.RDATA')
load('predABeveSIT_Y3.RDATA')
load('predABeveSIT_Y5.RDATA')
load('predABeveSIT_Y7.RDATA')
load('predABeveSIT_Y9.RDATA')
load('predABeveSIT_Y11.RDATA')
load('predABeveSIT_Y13.RDATA')
load('predABeveSIT_Y15.RDATA')
load('predABeveSIT_Y17.RDATA')
load('predABeveSIT_Y19.RDATA')
load('predABeveSIT_Y21.RDATA')
load('predABeveSIT_Y23.RDATA')
load('predABeveSIT_Y2.RDATA')
load('predABeveSIT_Y4.RDATA')
load('predABeveSIT_Y6.RDATA')
load('predABeveSIT_Y8.RDATA')
load('predABeveSIT_Y10.RDATA')
load('predABeveSIT_Y12.RDATA')
load('predABeveSIT_Y14.RDATA')
load('predABeveSIT_Y16.RDATA')
load('predABeveSIT_Y18.RDATA')
load('predABeveSIT_Y20.RDATA')
load('predABeveSIT_Y22.RDATA')


dat.iALLb = c(eveAll_Y1, eveAll_Y2, eveAll_Y3, eveAll_Y4, eveAll_Y5, eveAll_Y6, eveAll_Y7, eveAll_Y8, eveAll_Y9, eveAll_Y10, eveAll_Y11, 
              eveAll_Y12, eveAll_Y13, eveAll_Y14, eveAll_Y15, eveAll_Y16, eveAll_Y17, eveAll_Y18, eveAll_Y19, eveAll_Y20, eveAll_Y21, eveAll_Y22, eveAll_Y23)
dat.iALL = data.frame(eveAll = dat.iALLb, 
                      Year = c(rep(1998, length(eveAll_Y1)), rep(1999, length(eveAll_Y2)), rep(2000, length(eveAll_Y3)), rep(2001, length(eveAll_Y4)),
                               rep(2002, length(eveAll_Y5)),  rep(2003, length(eveAll_Y6)),rep(2004, length(eveAll_Y7)), rep(2005, length(eveAll_Y8)),
                               rep(2006, length(eveAll_Y9)),  rep(2007, length(eveAll_Y10)),rep(2008, length(eveAll_Y11)), rep(2009, length(eveAll_Y12)),
                               rep(2010, length(eveAll_Y13)),  rep(2011, length(eveAll_Y14)),rep(2012, length(eveAll_Y15)), rep(2013, length(eveAll_Y16)),
                               rep(2014, length(eveAll_Y17)),  rep(2015, length(eveAll_Y18)),rep(2016, length(eveAll_Y19)), rep(2017, length(eveAll_Y20)),
                               rep(2018, length(eveAll_Y21)),  rep(2019, length(eveAll_Y22)),rep(2020, length(eveAll_Y23))))

#####

#####
# Trend and residuals
######........................................................................................

##------------------ WL - traits

dat.iW1 = c(eveW1_Y1, eveW1_Y2, eveW1_Y3, eveW1_Y4, eveW1_Y5, eveW1_Y6, eveW1_Y7, eveW1_Y8, eveW1_Y9, eveW1_Y10, eveW1_Y11, 
            eveW1_Y12, eveW1_Y13, eveW1_Y14, eveW1_Y15, eveW1_Y16, eveW1_Y17, eveW1_Y18, eveW1_Y19, eveW1_Y20, eveW1_Y21, eveW1_Y22, eveW1_Y23)
dat.iW1 = data.frame(eveW1 = dat.iW1, 
                     Year = c(rep(1998, length(eveW1_Y1)), rep(1999, length(eveW1_Y2)), rep(2000, length(eveW1_Y3)), rep(2001, length(eveW1_Y4)),
                              rep(2002, length(eveW1_Y5)),  rep(2003, length(eveW1_Y6)),rep(2004, length(eveW1_Y7)), rep(2005, length(eveW1_Y8)),
                              rep(2006, length(eveW1_Y9)),  rep(2007, length(eveW1_Y10)),rep(2008, length(eveW1_Y11)), rep(2009, length(eveW1_Y12)),
                              rep(2010, length(eveW1_Y13)),  rep(2011, length(eveW1_Y14)),rep(2012, length(eveW1_Y15)), rep(2013, length(eveW1_Y16)),
                              rep(2014, length(eveW1_Y17)),  rep(2015, length(eveW1_Y18)),rep(2016, length(eveW1_Y19)), rep(2017, length(eveW1_Y20)),
                              rep(2018, length(eveW1_Y21)),  rep(2019, length(eveW1_Y22)),rep(2020, length(eveW1_Y23))))


dat.iW2 = c(eveW2_Y1, eveW2_Y2, eveW2_Y3, eveW2_Y4, eveW2_Y5, eveW2_Y6, eveW2_Y7, eveW2_Y8, eveW2_Y9, eveW2_Y10, eveW2_Y11, 
            eveW2_Y12, eveW2_Y13, eveW2_Y14, eveW2_Y15, eveW2_Y16, eveW2_Y17, eveW2_Y18, eveW2_Y19, eveW2_Y20, eveW2_Y21, eveW2_Y22, eveW2_Y23)
dat.iW2 = data.frame(eveW2 = dat.iW2, 
                     Year = c(rep(1998, length(eveW2_Y1)), rep(1999, length(eveW2_Y2)), rep(2000, length(eveW2_Y3)), rep(2001, length(eveW2_Y4)),
                              rep(2002, length(eveW2_Y5)),  rep(2003, length(eveW2_Y6)),rep(2004, length(eveW2_Y7)), rep(2005, length(eveW2_Y8)),
                              rep(2006, length(eveW2_Y9)),  rep(2007, length(eveW2_Y10)),rep(2008, length(eveW2_Y11)), rep(2009, length(eveW2_Y12)),
                              rep(2010, length(eveW2_Y13)),  rep(2011, length(eveW2_Y14)),rep(2012, length(eveW2_Y15)), rep(2013, length(eveW2_Y16)),
                              rep(2014, length(eveW2_Y17)),  rep(2015, length(eveW2_Y18)),rep(2016, length(eveW2_Y19)), rep(2017, length(eveW2_Y20)),
                              rep(2018, length(eveW2_Y21)),  rep(2019, length(eveW2_Y22)),rep(2020, length(eveW2_Y23))))


dat.iW3 = c(eveW3_Y1, eveW3_Y2, eveW3_Y3, eveW3_Y4, eveW3_Y5, eveW3_Y6, eveW3_Y7, eveW3_Y8, eveW3_Y9, eveW3_Y10, eveW3_Y11, 
            eveW3_Y12, eveW3_Y13, eveW3_Y14, eveW3_Y15, eveW3_Y16, eveW3_Y17, eveW3_Y18, eveW3_Y19, eveW3_Y20, eveW3_Y21, eveW3_Y22, eveW3_Y23)
dat.iW3 = data.frame(eveW3 = dat.iW3, 
                     Year = c(rep(1998, length(eveW3_Y1)), rep(1999, length(eveW3_Y2)), rep(2000, length(eveW3_Y3)), rep(2001, length(eveW3_Y4)),
                              rep(2002, length(eveW3_Y5)),  rep(2003, length(eveW3_Y6)),rep(2004, length(eveW3_Y7)), rep(2005, length(eveW3_Y8)),
                              rep(2006, length(eveW3_Y9)),  rep(2007, length(eveW3_Y10)),rep(2008, length(eveW3_Y11)), rep(2009, length(eveW3_Y12)),
                              rep(2010, length(eveW3_Y13)),  rep(2011, length(eveW3_Y14)),rep(2012, length(eveW3_Y15)), rep(2013, length(eveW3_Y16)),
                              rep(2014, length(eveW3_Y17)),  rep(2015, length(eveW3_Y18)),rep(2016, length(eveW3_Y19)), rep(2017, length(eveW3_Y20)),
                              rep(2018, length(eveW3_Y21)),  rep(2019, length(eveW3_Y22)),rep(2020, length(eveW3_Y23))))



ABeve = data.frame(eveH1=c(dat.iW1$eveW1), eveH2=c(dat.iW2$eveW2), eveH4=c(dat.iW3$eveW3),
                   eveall=c(dat.iALL$eveAll), SiteID = rep(Moth_knn2$sample, 23),
                   Year = c(rep(1998, length(eveW3_Y1)), rep(1999, length(eveW3_Y2)), rep(2000, length(eveW3_Y3)), rep(2001, length(eveW3_Y4)),
                            rep(2002, length(eveW3_Y5)),  rep(2003, length(eveW3_Y6)),rep(2004, length(eveW3_Y7)), rep(2005, length(eveW3_Y8)),
                            rep(2006, length(eveW3_Y9)),  rep(2007, length(eveW3_Y10)),rep(2008, length(eveW3_Y11)), rep(2009, length(eveW3_Y12)),
                            rep(2010, length(eveW3_Y13)),  rep(2011, length(eveW3_Y14)),rep(2012, length(eveW3_Y15)), rep(2013, length(eveW3_Y16)),
                            rep(2014, length(eveW3_Y17)),  rep(2015, length(eveW3_Y18)),rep(2016, length(eveW3_Y19)), rep(2017, length(eveW3_Y20)),
                            rep(2018, length(eveW3_Y21)),  rep(2019, length(eveW3_Y22)),rep(2020, length(eveW3_Y23))))


summary(ABeve)
head(ABeve)

#ABeve$Year = as.factor(ABeve$Year)
ABeve$Year = as.integer(ABeve$Year)
Site_un = Moth_knn2$sample


Trend_site_SITE = Resid_site_SITE = Resid2_site_SITE = coef_trend_SITE = matrix(NA, nrow= 109, ncol=4)
for (i in 1:109) {
  MapdataALL_i = ABeve[ABeve$SiteID == Site_un[i],]
  
  
  for (j in colnames(MapdataALL_i)[1:4]) {
    index_j = match(j, colnames(MapdataALL_i)[1:4])
    Fit_site_SITE = lm(reformulate('Year', j), data = MapdataALL_i)
    
    Resid_site_SITE[i, index_j] = var(resid(Fit_site_SITE))
    Resid2_site_SITE[i, index_j] = var(rstandard(Fit_site_SITE), na.rm=TRUE)
    coef_trend_SITE[i, index_j] = Fit_site_SITE$coefficients[2]
  }
}

Resid_site_SITE = as.data.frame(Resid_site_SITE)
colnames(Resid_site_SITE) = colnames(ABeve)[1:4]
Resid_site_SITE$spatial = 'SITE'

Resid2_site_SITE = as.data.frame(Resid2_site_SITE)
colnames(Resid2_site_SITE) = colnames(ABeve)[1:4]
Resid2_site_SITE$spatial = 'SITE'

coef_trend_SITE = as.data.frame(coef_trend_SITE)
colnames(coef_trend_SITE) = colnames(ABeve)[1:4]
coef_trend_SITE$spatial = 'SITE'

coef_trend_SITE$X = Moth_knn2$XEUREF
coef_trend_SITE$Y = Moth_knn2$YEUREF

Resid_site_SITE$X = Moth_knn2$XEUREF
Resid_site_SITE$Y = Moth_knn2$YEUREF

Resid2_site_SITE$X = Moth_knn2$XEUREF
Resid2_site_SITE$Y = Moth_knn2$YEUREF

## 
ABeve$X = Moth_knn2$XEUREF# same site order
ABeve$Y = Moth_knn2$YEUREF

## 

save(ABeve, Resid_site_SITE, Resid2_site_SITE, coef_trend_SITE,
     file='ABtrendWL_SITE_2808.RDATA')



##------------------ HP - traits

dat.iH1 = c(eveH1_Y1, eveH1_Y2, eveH1_Y3, eveH1_Y4, eveH1_Y5, eveH1_Y6, eveH1_Y7, eveH1_Y8, eveH1_Y9, eveH1_Y10, eveH1_Y11, 
            eveH1_Y12, eveH1_Y13, eveH1_Y14, eveH1_Y15, eveH1_Y16, eveH1_Y17, eveH1_Y18, eveH1_Y19, eveH1_Y20, eveH1_Y21, eveH1_Y22, eveH1_Y23)
dat.iH1 = data.frame(eveH1 = dat.iH1, 
                     Year = c(rep(1998, length(eveH1_Y1)), rep(1999, length(eveH1_Y2)), rep(2000, length(eveH1_Y3)), rep(2001, length(eveH1_Y4)),
                              rep(2002, length(eveH1_Y5)),  rep(2003, length(eveH1_Y6)),rep(2004, length(eveH1_Y7)), rep(2005, length(eveH1_Y8)),
                              rep(2006, length(eveH1_Y9)),  rep(2007, length(eveH1_Y10)),rep(2008, length(eveH1_Y11)), rep(2009, length(eveH1_Y12)),
                              rep(2010, length(eveH1_Y13)),  rep(2011, length(eveH1_Y14)),rep(2012, length(eveH1_Y15)), rep(2013, length(eveH1_Y16)),
                              rep(2014, length(eveH1_Y17)),  rep(2015, length(eveH1_Y18)),rep(2016, length(eveH1_Y19)), rep(2017, length(eveH1_Y20)),
                              rep(2018, length(eveH1_Y21)),  rep(2019, length(eveH1_Y22)),rep(2020, length(eveH1_Y23))))


dat.iH2 = c(eveH2_Y1, eveH2_Y2, eveH2_Y3, eveH2_Y4, eveH2_Y5, eveH2_Y6, eveH2_Y7, eveH2_Y8, eveH2_Y9, eveH2_Y10, eveH2_Y11, 
            eveH2_Y12, eveH2_Y13, eveH2_Y14, eveH2_Y15, eveH2_Y16, eveH2_Y17, eveH2_Y18, eveH2_Y19, eveH2_Y20, eveH2_Y21, eveH2_Y22, eveH2_Y23)
dat.iH2 = data.frame(eveH2 = dat.iH2, 
                     Year = c(rep(1998, length(eveH2_Y1)), rep(1999, length(eveH2_Y2)), rep(2000, length(eveH2_Y3)), rep(2001, length(eveH2_Y4)),
                              rep(2002, length(eveH2_Y5)),  rep(2003, length(eveH2_Y6)),rep(2004, length(eveH2_Y7)), rep(2005, length(eveH2_Y8)),
                              rep(2006, length(eveH2_Y9)),  rep(2007, length(eveH2_Y10)),rep(2008, length(eveH2_Y11)), rep(2009, length(eveH2_Y12)),
                              rep(2010, length(eveH2_Y13)),  rep(2011, length(eveH2_Y14)),rep(2012, length(eveH2_Y15)), rep(2013, length(eveH2_Y16)),
                              rep(2014, length(eveH2_Y17)),  rep(2015, length(eveH2_Y18)),rep(2016, length(eveH2_Y19)), rep(2017, length(eveH2_Y20)),
                              rep(2018, length(eveH2_Y21)),  rep(2019, length(eveH2_Y22)),rep(2020, length(eveH2_Y23))))


dat.iH4 = c(eveH4_Y1, eveH4_Y2, eveH4_Y3, eveH4_Y4, eveH4_Y5, eveH4_Y6, eveH4_Y7, eveH4_Y8, eveH4_Y9, eveH4_Y10, eveH4_Y11, 
            eveH4_Y12, eveH4_Y13, eveH4_Y14, eveH4_Y15, eveH4_Y16, eveH4_Y17, eveH4_Y18, eveH4_Y19, eveH4_Y20, eveH4_Y21, eveH4_Y22, eveH4_Y23)
dat.iH4 = data.frame(eveH4 = dat.iH4, 
                     Year = c(rep(1998, length(eveH4_Y1)), rep(1999, length(eveH4_Y2)), rep(2000, length(eveH4_Y3)), rep(2001, length(eveH4_Y4)),
                              rep(2002, length(eveH4_Y5)),  rep(2003, length(eveH4_Y6)),rep(2004, length(eveH4_Y7)), rep(2005, length(eveH4_Y8)),
                              rep(2006, length(eveH4_Y9)),  rep(2007, length(eveH4_Y10)),rep(2008, length(eveH4_Y11)), rep(2009, length(eveH4_Y12)),
                              rep(2010, length(eveH4_Y13)),  rep(2011, length(eveH4_Y14)),rep(2012, length(eveH4_Y15)), rep(2013, length(eveH4_Y16)),
                              rep(2014, length(eveH4_Y17)),  rep(2015, length(eveH4_Y18)),rep(2016, length(eveH4_Y19)), rep(2017, length(eveH4_Y20)),
                              rep(2018, length(eveH4_Y21)),  rep(2019, length(eveH4_Y22)),rep(2020, length(eveH4_Y23))))


dat.iH31 = c(eveH31_Y1, eveH31_Y2, eveH31_Y3, eveH31_Y4, eveH31_Y5, eveH31_Y6, eveH31_Y7, eveH31_Y8, eveH31_Y9, eveH31_Y10, eveH31_Y11, 
             eveH31_Y12, eveH31_Y13, eveH31_Y14, eveH31_Y15, eveH31_Y16, eveH31_Y17, eveH31_Y18, eveH31_Y19, eveH31_Y20, eveH31_Y21, eveH31_Y22, eveH31_Y23)
dat.iH31 = data.frame(eveH31 = dat.iH31, 
                      Year = c(rep(1998, length(eveH31_Y1)), rep(1999, length(eveH31_Y2)), rep(2000, length(eveH31_Y3)), rep(2001, length(eveH31_Y4)),
                               rep(2002, length(eveH31_Y5)),  rep(2003, length(eveH31_Y6)),rep(2004, length(eveH31_Y7)), rep(2005, length(eveH31_Y8)),
                               rep(2006, length(eveH31_Y9)),  rep(2007, length(eveH31_Y10)),rep(2008, length(eveH31_Y11)), rep(2009, length(eveH31_Y12)),
                               rep(2010, length(eveH31_Y13)),  rep(2011, length(eveH31_Y14)),rep(2012, length(eveH31_Y15)), rep(2013, length(eveH31_Y16)),
                               rep(2014, length(eveH31_Y17)),  rep(2015, length(eveH31_Y18)),rep(2016, length(eveH31_Y19)), rep(2017, length(eveH31_Y20)),
                               rep(2018, length(eveH31_Y21)),  rep(2019, length(eveH31_Y22)),rep(2020, length(eveH31_Y23))))



dat.iH32 = c(eveH32_Y1, eveH32_Y2, eveH32_Y3, eveH32_Y4, eveH32_Y5, eveH32_Y6, eveH32_Y7, eveH32_Y8, eveH32_Y9, eveH32_Y10, eveH32_Y11, 
             eveH32_Y12, eveH32_Y13, eveH32_Y14, eveH32_Y15, eveH32_Y16, eveH32_Y17, eveH32_Y18, eveH32_Y19, eveH32_Y20, eveH32_Y21, eveH32_Y22, eveH32_Y23)
dat.iH32 = data.frame(eveH32 = dat.iH32, 
                      Year = c(rep(1998, length(eveH32_Y1)), rep(1999, length(eveH32_Y2)), rep(2000, length(eveH32_Y3)), rep(2001, length(eveH32_Y4)),
                               rep(2002, length(eveH32_Y5)),  rep(2003, length(eveH32_Y6)),rep(2004, length(eveH32_Y7)), rep(2005, length(eveH32_Y8)),
                               rep(2006, length(eveH32_Y9)),  rep(2007, length(eveH32_Y10)),rep(2008, length(eveH32_Y11)), rep(2009, length(eveH32_Y12)),
                               rep(2010, length(eveH32_Y13)),  rep(2011, length(eveH32_Y14)),rep(2012, length(eveH32_Y15)), rep(2013, length(eveH32_Y16)),
                               rep(2014, length(eveH32_Y17)),  rep(2015, length(eveH32_Y18)),rep(2016, length(eveH32_Y19)), rep(2017, length(eveH32_Y20)),
                               rep(2018, length(eveH32_Y21)),  rep(2019, length(eveH32_Y22)),rep(2020, length(eveH32_Y23))))


dat.iH33 = c(eveH33_Y1, eveH33_Y2, eveH33_Y3, eveH33_Y4, eveH33_Y5, eveH33_Y6, eveH33_Y7, eveH33_Y8, eveH33_Y9, eveH33_Y10, eveH33_Y11, 
             eveH33_Y12, eveH33_Y13, eveH33_Y14, eveH33_Y15, eveH33_Y16, eveH33_Y17, eveH33_Y18, eveH33_Y19, eveH33_Y20, eveH33_Y21, eveH33_Y22, eveH33_Y23)
dat.iH33 = data.frame(eveH33 = dat.iH33, 
                      Year = c(rep(1998, length(eveH33_Y1)), rep(1999, length(eveH33_Y2)), rep(2000, length(eveH33_Y3)), rep(2001, length(eveH33_Y4)),
                               rep(2002, length(eveH33_Y5)),  rep(2003, length(eveH33_Y6)),rep(2004, length(eveH33_Y7)), rep(2005, length(eveH33_Y8)),
                               rep(2006, length(eveH33_Y9)),  rep(2007, length(eveH33_Y10)),rep(2008, length(eveH33_Y11)), rep(2009, length(eveH33_Y12)),
                               rep(2010, length(eveH33_Y13)),  rep(2011, length(eveH33_Y14)),rep(2012, length(eveH33_Y15)), rep(2013, length(eveH33_Y16)),
                               rep(2014, length(eveH33_Y17)),  rep(2015, length(eveH33_Y18)),rep(2016, length(eveH33_Y19)), rep(2017, length(eveH33_Y20)),
                               rep(2018, length(eveH33_Y21)),  rep(2019, length(eveH33_Y22)),rep(2020, length(eveH33_Y23))))


ABeve = data.frame(eveH1=c(dat.iH1$eveH1), eveH2=c(dat.iH2$eveH2), eveH4=c(dat.iH4$eveH4),
                   eveH31=c(dat.iH31$eveH31), eveH32=c(dat.iH32$eveH32), eveH33=c(dat.iH33$eveH33),
                   eveall=c(dat.iALL$eveAll), SiteID = rep(Moth_knn2$sample, 23),
                   Year = c(rep(1998, length(eveH33_Y1)), rep(1999, length(eveH33_Y2)), rep(2000, length(eveH33_Y3)), rep(2001, length(eveH33_Y4)),
                            rep(2002, length(eveH33_Y5)),  rep(2003, length(eveH33_Y6)),rep(2004, length(eveH33_Y7)), rep(2005, length(eveH33_Y8)),
                            rep(2006, length(eveH33_Y9)),  rep(2007, length(eveH33_Y10)),rep(2008, length(eveH33_Y11)), rep(2009, length(eveH33_Y12)),
                            rep(2010, length(eveH33_Y13)),  rep(2011, length(eveH33_Y14)),rep(2012, length(eveH33_Y15)), rep(2013, length(eveH33_Y16)),
                            rep(2014, length(eveH33_Y17)),  rep(2015, length(eveH33_Y18)),rep(2016, length(eveH33_Y19)), rep(2017, length(eveH33_Y20)),
                            rep(2018, length(eveH33_Y21)),  rep(2019, length(eveH33_Y22)),rep(2020, length(eveH33_Y23))))


summary(ABeve)
head(ABeve)
ABeve$Year = as.integer(ABeve$Year)
Site_un = Moth_knn2$sample


Trend_site_SITE = Resid_site_SITE = Resid2_site_SITE = coef_trend_SITE = matrix(NA, nrow= 109, ncol=7)
for (i in 1:109) {
  MapdataALL_i = ABeve[ABeve$SiteID == Site_un[i],]
  
  
  for (j in colnames(MapdataALL_i)[1:7]) {
    index_j = match(j, colnames(MapdataALL_i)[1:7])
    Fit_site_SITE = lm(reformulate('Year', j), data = MapdataALL_i)
    
    Resid_site_SITE[i, index_j] = var(resid(Fit_site_SITE))
    Resid2_site_SITE[i, index_j] = var(rstandard(Fit_site_SITE), na.rm=TRUE)
    coef_trend_SITE[i, index_j] = Fit_site_SITE$coefficients[2]
  }
}

Resid_site_SITE = as.data.frame(Resid_site_SITE)
colnames(Resid_site_SITE) = colnames(ABeve)[1:7]
Resid_site_SITE$spatial = 'SITE'

Resid2_site_SITE = as.data.frame(Resid2_site_SITE)
colnames(Resid2_site_SITE) = colnames(ABeve)[1:7]
Resid2_site_SITE$spatial = 'SITE'

coef_trend_SITE = as.data.frame(coef_trend_SITE)
colnames(coef_trend_SITE) = colnames(ABeve)[1:7]
coef_trend_SITE$spatial = 'SITE'


coef_trend_SITE$X = Moth_knn2$XEUREF
coef_trend_SITE$Y = Moth_knn2$YEUREF

Resid_site_SITE$X = Moth_knn2$XEUREF
Resid_site_SITE$Y = Moth_knn2$YEUREF

Resid2_site_SITE$X = Moth_knn2$XEUREF
Resid2_site_SITE$Y = Moth_knn2$YEUREF

## 
ABeve$X = Moth_knn2$XEUREF# same site order
ABeve$Y = Moth_knn2$YEUREF

## 

save(ABeve, Resid_site_SITE, Resid2_site_SITE, coef_trend_SITE,
     file='ABtrendHP_SITE_2808.RDATA')
