##################################################################################
#
#		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
#........................................................................................
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/1.EnvironmentalInfo")
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('predPArichSIT_Y1.RDATA')
load('predPArichSIT_Y3.RDATA')
load('predPArichSIT_Y5.RDATA')
load('predPArichSIT_Y7.RDATA')
load('predPArichSIT_Y9.RDATA')
load('predPArichSIT_Y11.RDATA')
load('predPArichSIT_Y13.RDATA')
load('predPArichSIT_Y15.RDATA')
load('predPArichSIT_Y17.RDATA')
load('predPArichSIT_Y19.RDATA')
load('predPArichSIT_Y21.RDATA')
load('predPArichSIT_Y23.RDATA')

load('predPArichSIT_Y2.RDATA')
load('predPArichSIT_Y4.RDATA')
load('predPArichSIT_Y6.RDATA')
load('predPArichSIT_Y8.RDATA')
load('predPArichSIT_Y10.RDATA')
load('predPArichSIT_Y12.RDATA')
load('predPArichSIT_Y14.RDATA')
load('predPArichSIT_Y16.RDATA')
load('predPArichSIT_Y18.RDATA')
load('predPArichSIT_Y20.RDATA')
load('predPArichSIT_Y22.RDATA')


## All

dat.iALLb = c(S_Y1, S_Y2, S_Y3, S_Y4, S_Y5, S_Y6, S_Y7, S_Y8, S_Y9, S_Y10, S_Y11, 
              S_Y12, S_Y13, S_Y14, S_Y15, S_Y16, S_Y17, S_Y18, S_Y19, S_Y20, S_Y21, S_Y22, S_Y23)
dat.iALL = data.frame(S = dat.iALLb, 
                      Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))




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

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

dat.iW1 = c(Raw1HW1_Y1, Raw1HW1_Y2, Raw1HW1_Y3, Raw1HW1_Y4, Raw1HW1_Y5, Raw1HW1_Y6, Raw1HW1_Y7, Raw1HW1_Y8, Raw1HW1_Y9, Raw1HW1_Y10, Raw1HW1_Y11, 
            Raw1HW1_Y12, Raw1HW1_Y13, Raw1HW1_Y14, Raw1HW1_Y15, Raw1HW1_Y16, Raw1HW1_Y17, Raw1HW1_Y18, Raw1HW1_Y19, Raw1HW1_Y20, Raw1HW1_Y21, Raw1HW1_Y22, Raw1HW1_Y23)
dat.iW1 = data.frame(S = dat.iW1, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                              rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                              rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                              rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                              rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                              rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



dat.iW2 = c(Raw1HW2_Y1, Raw1HW2_Y2, Raw1HW2_Y3, Raw1HW2_Y4, Raw1HW2_Y5, Raw1HW2_Y6, Raw1HW2_Y7, Raw1HW2_Y8, Raw1HW2_Y9, Raw1HW2_Y10, Raw1HW2_Y11, 
            Raw1HW2_Y12, Raw1HW2_Y13, Raw1HW2_Y14, Raw1HW2_Y15, Raw1HW2_Y16, Raw1HW2_Y17, Raw1HW2_Y18, Raw1HW2_Y19, Raw1HW2_Y20, Raw1HW2_Y21, Raw1HW2_Y22, Raw1HW2_Y23)
dat.iW2 = data.frame(S = dat.iW2, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                              rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                              rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                              rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                              rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                              rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))

dat.iW3 = c(Raw1HW3_Y1, Raw1HW3_Y2, Raw1HW3_Y3, Raw1HW3_Y4, Raw1HW3_Y5, Raw1HW3_Y6, Raw1HW3_Y7, Raw1HW3_Y8, Raw1HW3_Y9, Raw1HW3_Y10, Raw1HW3_Y11, 
            Raw1HW3_Y12, Raw1HW3_Y13, Raw1HW3_Y14, Raw1HW3_Y15, Raw1HW3_Y16, Raw1HW3_Y17, Raw1HW3_Y18, Raw1HW3_Y19, Raw1HW3_Y20, Raw1HW3_Y21, Raw1HW3_Y22, Raw1HW3_Y23)
dat.iW3 = data.frame(S = dat.iW3, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                              rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                              rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                              rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                              rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                              rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



PArich = data.frame(RicH1=c(dat.iW1$S), RicH2=c(dat.iW2$S), RicH4=c(dat.iW3$S),
                    Richall=c(dat.iALL$S), SiteID = rep(Moth_knn2$sample, 23),
                    Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                             rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                             rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                             rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                             rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                             rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))


summary(PArich)
head(PArich)

Site_un = Moth_knn2$sample
PArich$Year = as.integer(PArich$Year)

Trend_site_SITE = Resid_site_SITE = Resid2_site_SITE = coef_trend_SITE = matrix(NA, nrow= 109, ncol=4)
for (i in 1:109) {
  PArich_i = PArich[PArich$SiteID == Site_un[i],]
  
  for (j in colnames(PArich_i)[1:4]) {
    index_j = match(j, colnames(PArich_i)[1:4])
    Fit_site_SITE = lm(reformulate('Year', j), data = PArich_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(PArich)[1:4]
Resid_site_SITE$spatial = 'SITE'

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

coef_trend_SITE = as.data.frame(coef_trend_SITE)
colnames(coef_trend_SITE) = colnames(PArich)[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

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

save(PArich, Resid_site_SITE, Resid2_site_SITE, coef_trend_SITE,
     file='PAtrendWL_SITE_2808.RDATA')



##------------------ HP - trait

dat.iH1 = c(Raw1H1_Y1, Raw1H1_Y2, Raw1H1_Y3, Raw1H1_Y4, Raw1H1_Y5, Raw1H1_Y6, Raw1H1_Y7, Raw1H1_Y8, Raw1H1_Y9, Raw1H1_Y10, Raw1H1_Y11, 
            Raw1H1_Y12, Raw1H1_Y13, Raw1H1_Y14, Raw1H1_Y15, Raw1H1_Y16, Raw1H1_Y17, Raw1H1_Y18, Raw1H1_Y19, Raw1H1_Y20, Raw1H1_Y21, Raw1H1_Y22, Raw1H1_Y23)
dat.iH1 = data.frame(S = dat.iH1, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                              rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                              rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                              rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                              rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                              rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



dat.iH2 = c(Raw1H2_Y1, Raw1H2_Y2, Raw1H2_Y3, Raw1H2_Y4, Raw1H2_Y5, Raw1H2_Y6, Raw1H2_Y7, Raw1H2_Y8, Raw1H2_Y9, Raw1H2_Y10, Raw1H2_Y11, 
            Raw1H2_Y12, Raw1H2_Y13, Raw1H2_Y14, Raw1H2_Y15, Raw1H2_Y16, Raw1H2_Y17, Raw1H2_Y18, Raw1H2_Y19, Raw1H2_Y20, Raw1H2_Y21, Raw1H2_Y22, Raw1H2_Y23)
dat.iH2 = data.frame(S = dat.iH2, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                              rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                              rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                              rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                              rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                              rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))


dat.iH4 = c(Raw1H4_Y1, Raw1H4_Y2, Raw1H4_Y3, Raw1H4_Y4, Raw1H4_Y5, Raw1H4_Y6, Raw1H4_Y7, Raw1H4_Y8, Raw1H4_Y9, Raw1H4_Y10, Raw1H4_Y11, 
            Raw1H4_Y12, Raw1H4_Y13, Raw1H4_Y14, Raw1H4_Y15, Raw1H4_Y16, Raw1H4_Y17, Raw1H4_Y18, Raw1H4_Y19, Raw1H4_Y20, Raw1H4_Y21, Raw1H4_Y22, Raw1H4_Y23)
dat.iH4 = data.frame(S = dat.iH4, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                              rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                              rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                              rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                              rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                              rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



dat.iH31 = c(Raw1H31_Y1, Raw1H31_Y2, Raw1H31_Y3, Raw1H31_Y4, Raw1H31_Y5, Raw1H31_Y6, Raw1H31_Y7, Raw1H31_Y8, Raw1H31_Y9, Raw1H31_Y10, Raw1H31_Y11, 
             Raw1H31_Y12, Raw1H31_Y13, Raw1H31_Y14, Raw1H31_Y15, Raw1H31_Y16, Raw1H31_Y17, Raw1H31_Y18, Raw1H31_Y19, Raw1H31_Y20, Raw1H31_Y21, Raw1H31_Y22, Raw1H31_Y23)
dat.iH31 = data.frame(S = dat.iH31, 
                      Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



dat.iH32 = c(Raw1H32_Y1, Raw1H32_Y2, Raw1H32_Y3, Raw1H32_Y4, Raw1H32_Y5, Raw1H32_Y6, Raw1H32_Y7, Raw1H32_Y8, Raw1H32_Y9, Raw1H32_Y10, Raw1H32_Y11, 
             Raw1H32_Y12, Raw1H32_Y13, Raw1H32_Y14, Raw1H32_Y15, Raw1H32_Y16, Raw1H32_Y17, Raw1H32_Y18, Raw1H32_Y19, Raw1H32_Y20, Raw1H32_Y21, Raw1H32_Y22, Raw1H32_Y23)
dat.iH32 = data.frame(S = dat.iH32, 
                      Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



dat.iH33 = c(Raw1H33_Y1, Raw1H33_Y2, Raw1H33_Y3, Raw1H33_Y4, Raw1H33_Y5, Raw1H33_Y6, Raw1H33_Y7, Raw1H33_Y8, Raw1H33_Y9, Raw1H33_Y10, Raw1H33_Y11, 
             Raw1H33_Y12, Raw1H33_Y13, Raw1H33_Y14, Raw1H33_Y15, Raw1H33_Y16, Raw1H33_Y17, Raw1H33_Y18, Raw1H33_Y19, Raw1H33_Y20, Raw1H33_Y21, Raw1H33_Y22, Raw1H33_Y23)
dat.iH33 = data.frame(S = dat.iH33, 
                      Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



PArich = data.frame(RicH1=c(dat.iH1$S), RicH2=c(dat.iH2$S), RicH4=c(dat.iH4$S),
                    RicH31=c(dat.iH31$S), RicH32=c(dat.iH32$S), RicH33=c(dat.iH33$S),
                    Richall=c(dat.iALL$S), SiteID = rep(Moth_knn2$sample, 23),
                    Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                             rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                             rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                             rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                             rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                             rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))

summary(PArich)
head(PArich)

Site_un = Moth_knn2$sample
PArich$Year = as.integer(PArich$Year)

Trend_site_SITE = Resid_site_SITE = Resid2_site_SITE = coef_trend_SITE = matrix(NA, nrow= 109, ncol=7)
for (i in 1:109) {
  PArich_i = PArich[PArich$SiteID == Site_un[i],]
  
  for (j in colnames(PArich_i)[1:7]) {
    index_j = match(j, colnames(PArich_i)[1:7])
    Fit_site_SITE = lm(reformulate('Year', j), data = PArich_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(PArich)[1:7]
Resid_site_SITE$spatial = 'SITE'

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

coef_trend_SITE = as.data.frame(coef_trend_SITE)
colnames(coef_trend_SITE) = colnames(PArich)[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

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

save(PArich, Resid_site_SITE, Resid2_site_SITE, coef_trend_SITE,
     file='PAtrendHP_SITE_2808.RDATA')

