mod8p3

bac<-read.table("bac_dust.txt", header = T, sep = "\t")
library(maps)

## Warning: package 'maps' was built under R version 4.3.2

library(ggplot2)
world_map <- map_data("world")
library(ggplot2)
ggplot(world_map)+
  geom_polygon(aes(x = long, y = lat, group = group), fill = "lightgray", colour = "black", size = 0.1)+
  geom_point(data = bac, aes(x = Longitude, y = Latitude, color = Continent),  size = 2)+
  scale_color_manual(values = c("SouthAmerica" = 'chartreuse4', "NorthAmerica" = 'chartreuse1', "Africa" = 'goldenrod1', "Oceania" = 'darkorchid3', "Europe" = 'firebrick4', "Asia" = 'cornflowerblue'))+
  theme_bw(base_size = 15)+
  theme(legend.title = element_blank(), axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank())

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ggplot(world_map)+
  geom_polygon(aes(x = long, y = lat, group = group), fill = "lightgray", colour = "black", size = 0.1)+
  geom_point(data = bac, aes(x = Longitude, y = Latitude, color = Continent, size=Richness), alpha=0.6)+
  scale_color_manual(values = c("SouthAmerica" = 'chartreuse4', "NorthAmerica" = 'chartreuse1', "Africa" = 'goldenrod1', "Oceania" = 'darkorchid3', "Europe" = 'firebrick4', "Asia" = 'cornflowerblue'))+
  theme_bw(base_size = 15)+
  theme(legend.title = element_blank(), axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank())

bac.dists <- as.matrix(dist(cbind(bac$Latitude, bac$Longitude)))
bac.dists.inv <- 1/bac.dists
diag(bac.dists.inv) <- 0
library(ape)

## Warning: package 'ape' was built under R version 4.3.2

Moran.I(bac$Richness, bac.dists.inv)

## $observed
## [1] 0.1159867
## 
## $expected
## [1] -0.003225806
## 
## $sd
## [1] 0.03860909
## 
## $p.value
## [1] 0.002017263

Moran.I(bac$Proteobacteria, bac.dists.inv)

## $observed
## [1] 0.216501
## 
## $expected
## [1] -0.003225806
## 
## $sd
## [1] 0.03858986
## 
## $p.value
## [1] 1.241692e-08

anova(lm(Richness~Proteobacteria, data=bac)) #Not significant

## Analysis of Variance Table
## 
## Response: Richness
##                 Df   Sum Sq Mean Sq F value Pr(>F)
## Proteobacteria   1    38901   38901  0.5349 0.4651
## Residuals      309 22471204   72722

bac$resid <-resid(lm(Richness~Proteobacteria, data=bac))

ggplot(world_map)+
  geom_polygon(aes(x = long, y = lat, group = group), fill = "lightgray", colour = "black", size = 0.1)+
  geom_point(data = bac, aes(x = Longitude, y = Latitude, color = Continent, size=resid), alpha=0.6)+
  scale_color_manual(values = c("SouthAmerica" = 'chartreuse4', "NorthAmerica" = 'chartreuse1', "Africa" = 'goldenrod1', "Oceania" = 'darkorchid3', "Europe" = 'firebrick4', "Asia" = 'cornflowerblue'))+
  theme_bw(base_size = 15)+
  theme(legend.title = element_blank(), axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank())

Moran.I(bac$resid, bac.dists.inv)

## $observed
## [1] 0.1098749
## 
## $expected
## [1] -0.003225806
## 
## $sd
## [1] 0.03861023
## 
## $p.value
## [1] 0.003397342

library(spaMM)

## Warning: package 'spaMM' was built under R version 4.3.2

## Registered S3 methods overwritten by 'registry':
##   method               from 
##   print.registry_field proxy
##   print.registry_entry proxy

## spaMM (Rousset & Ferdy, 2014, version 4.4.0) is loaded.
## Type 'help(spaMM)' for a short introduction,
## 'news(package='spaMM')' for news,
## and 'citation('spaMM')' for proper citation.
## Further infos, slides, etc. at https://gitlab.mbb.univ-montp2.fr/francois/spamm-ref.

fitme(Richness~Proteobacteria+Matern(1|Latitude+Longitude), data = bac)

## If the 'RSpectra' package were installed, an extreme eigenvalue computation could be faster.

## formula: Richness ~ Proteobacteria + Matern(1 | Latitude + Longitude)
## ML: Estimation of corrPars, lambda and phi by ML.
##     Estimation of fixed effects by ML.
## Estimation of lambda and phi by 'outer' ML, maximizing logL.
## family: gaussian( link = identity ) 
##  ------------ Fixed effects (beta) ------------
##                Estimate Cond. SE t-value
## (Intercept)      411.58    48.47  8.4910
## Proteobacteria    37.66   112.90  0.3336
##  --------------- Random effects ---------------
## Family: gaussian( link = identity ) 
##                    --- Correlation parameters:
##      1.nu     1.rho 
## 0.2041736 0.3145350 
##            --- Variance parameters ('lambda'):
## lambda = var(u) for u ~ Gaussian; 
##    Latitude .  :  19190  
## # of obs: 311; # of groups: Latitude ., 311 
##  -------------- Residual variance  ------------
## phi estimate was 54713 
##  ------------- Likelihood values  -------------
##                         logLik
## logL       (p_v(h)): -2170.262

#nu     rho 
#0.2041723 0.3145336
dd<-dist(bac[,c("Longitude", "Latitude")])
mm<-MaternCorr(dd, nu=0.2041723, rho=0.3145336)
plot(as.numeric(dd), as.numeric(mm), xlab="Distance between pairs of locations", ylab="Estimated correlation")