library(readxl)
locations <- read_excel("locations.xlsx")

To make maps on R, you need to register with Google Cloud Platform (free trial) and get an API key. Walkthrough for setting this up is here -> https://www.littlemissdata.com/blog/maps

Please note that the above walkthrough does not tell you to register Geocoding, only API Maps Static. You need both to do what I have done here. It will give you one key, with both registered to it.

You need to load ggmap and register your google API key for API Maps Static and Geocoding. You also need an internet connection to knit it as it is using google maps to find the coordinates.

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
ggmap::register_google(key = "AIzaSyB-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0") #your API for Maps Static and Geocoding on Google Cloud. It must be set to unrestricted in Google.
geo <- geocode("Sydney", size = 0.3) #find coordinates of where you want the centre of your map to be.
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Sydney&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
print(geo)
## # A tibble: 1 x 2
##     lon   lat
##   <dbl> <dbl>
## 1  151. -33.9
register_google(key = "AIzaSyB-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0")
sc <- mutate_geocode(locations, SC) #calculate coordinates for each location. syntax = mutate_geocode(#YourDatset, #Your Locations Variable Name) 
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Westfield+Eastgardens&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Fashion+Spree+Orange+Grove&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Neeta+City&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Marrickville+Metro&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Westfield+Penrith&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Westpoint+Shopping+Centre,+Blacktown&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Randwick+Shopping+Centre&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Kellyville+Village+Shops&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Rouse+Hill+Town+Centre&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Hurstville+Westfield&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Revesby+Village+Centre&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1-3+Brady+Street,+Mosman+NSW+2088&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Birkenhead+Point&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Westfield+Hornsby&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Minto+Marketplace&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Castle+Towers+Shopping+Centre,+6-14+Castle+Street,+Castle+Hill+NSW+2154&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Westfield+Parramatta&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Westfield+Bondi+Junction&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Top+Ryde+City+Shopping+Centre&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Rhodes+Waterside&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Westfield+Burwood&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=DFO+Homebush&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Bankstown+Centro&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Ashfield+Mall&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Macquarie+Centre&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Warringah+Mall&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Westield+Miranda&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Dee+Why+Grand&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## "Dee Why Grand" not uniquely geocoded, using "19 pacific parade, dee why nsw 2099, australia"
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Westfield+Chatswood&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Stockland+Piccadilly+Shopping+Centre&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Westfield+Mount+Druitt&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=South+Village,+580+Princes+Highway+Kirrawee+NSW+2232&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1+Bay+Street,+Ultimo+NSW+2007&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Strathfield+Plaza&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Kings+Cross+Centre&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Auburn+Cebtral&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ tibble  3.0.4     ✓ purrr   0.3.4
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
sc2 <- as_tibble(sc) #make geocoded data into its own dataset
print(sc2, n = 35) #check to see that it is correct. sometimes it will put down the address of a similar named place in another location. Check to see if any values are very different, and then change your dataset to be more specific (like a full street address)
## # A tibble: 36 x 9
##    SC             Suburb  Person Population Income   Pop IncomeValue   lon   lat
##    <chr>          <chr>   <chr>  <chr>      <chr>  <dbl>       <dbl> <dbl> <dbl>
##  1 Westfield Eas… Eastga… Tarryn LOW        HIGH     860        1833  151. -33.9
##  2 Fashion Spree… Liverp… Tarryn HIGH       LOW    27804        1089  151. -33.9
##  3 Neeta City     Fairfi… Tarryn LOW        LOW    18081         975  151. -33.9
##  4 Marrickville … Marric… Tarryn HIGH       HIGH   26592        1814  151. -33.9
##  5 Westfield Pen… Penrith Gemma  LOW        LOW    10468        2559  151. -33.8
##  6 Westpoint Sho… Blackt… Gemma  HIGH       LOW    47176        1459  151. -33.8
##  7 Randwick Shop… Randwi… Gemma  LOW        LOW    13295        1142  151. -33.9
##  8 Kellyville Vi… Kellyv… Gemma  HIGH       HIGH   29986        2169  151. -33.7
##  9 Rouse Hill To… Rouse … Charl… LOW        HIGH    7965        2401  151. -33.7
## 10 Hurstville We… Hurstv… Charl… HIGH       LOW    29822        1382  151. -34.0
## 11 Revesby Villa… Revesby Charl… LOW        LOW    14176        1516  151. -34.0
## 12 1-3 Brady Str… Mosman  Charl… HIGH       HIGH   28475        2533  151. -33.8
## 13 Birkenhead Po… Drummo… Corne… LOW        HIGH   11950        2353  151. -33.9
## 14 Westfield Hor… Hornsby Corne… HIGH       LOW    12551        1383  151. -33.7
## 15 Minto Marketp… Minto   Corne… LOW        LOW    22168        1696  151. -34.0
## 16 Castle Towers… Castle… Corne… HIGH       HIGH   39594        2219  151. -33.7
## 17 Westfield Par… Parram… Yejin  LOW        HIGH   25798        1739  151. -33.8
## 18 Westfield Bon… Bondi … Yejin  HIGH       LOW    10986        2026  151. -33.9
## 19 Top Ryde City… Top Ry… Yejin  LOW        LOW     7929        1604  151. -33.8
## 20 Rhodes Waters… Rhodes  Yejin  HIGH       HIGH   13232        1617  151. -33.8
## 21 Westfield Bur… Burwood Domin… LOW        LOW    16030        1398  151. -33.9
## 22 DFO Homebush   Homebu… Domin… LOW        HIGH    7007        1810  151. -33.9
## 23 Bankstown Cen… Bankst… Domin… HIGH       LOW    32113        1120  151. -33.9
## 24 Ashfield Mall  Ashfie… Domin… HIGH       HIGH   23841        1632  151. -33.9
## 25 Macquarie Cen… Macqua… Ella   LOW        LOW     8144        1631  151. -33.8
## 26 Warringah Mall Brookv… Ella   LOW        HIGH    3161        1822  151. -33.8
## 27 Westield Mira… Miranda Ella   HIGH       LOW    29653        1267  151. -34.0
## 28 Dee Why Grand  Dee Why Ella   HIGH       HIGH   21518        1716  151. -33.8
## 29 Westfield Cha… Chatsw… Kiki   HIGH       HIGH   24913        1881  151. -33.8
## 30 Stockland Pic… Sydney  Kiki   LOW        HIGH   17252        1949  151. -33.9
## 31 Westfield Mou… Mount … Kiki   LOW        LOW    16726        1268  151. -33.8
## 32 South Village… Suther… Kiki   HIGH       LOW    20868        1692  151. -34.0
## 33 1 Bay Street,… Ultimo  Vivvy  LOW        LOW     8845        1230  151. -33.9
## 34 Strathfield P… Strath… Vivvy  HIGH       HIGH   25813        1892  151. -33.9
## 35 Kings Cross C… Potts … Vivvy  LOW        HIGH   17252        1949  151. -33.9
## # … with 1 more row
str(sc2)
## tibble [36 × 9] (S3: tbl_df/tbl/data.frame)
##  $ SC         : chr [1:36] "Westfield Eastgardens" "Fashion Spree Orange Grove" "Neeta City" "Marrickville Metro" ...
##  $ Suburb     : chr [1:36] "Eastgardens" "Liverpool" "Fairfield" "Marrickville" ...
##  $ Person     : chr [1:36] "Tarryn" "Tarryn" "Tarryn" "Tarryn" ...
##  $ Population : chr [1:36] "LOW" "HIGH" "LOW" "HIGH" ...
##  $ Income     : chr [1:36] "HIGH" "LOW" "LOW" "HIGH" ...
##  $ Pop        : num [1:36] 860 27804 18081 26592 10468 ...
##  $ IncomeValue: num [1:36] 1833 1089 975 1814 2559 ...
##  $ lon        : num [1:36] 151 151 151 151 151 ...
##  $ lat        : num [1:36] -33.9 -33.9 -33.9 -33.9 -33.8 ...

mapview


library(ggplot2)





```r
geocode("Parramatta") #check coordinates of where you want your centre of the map to be. I changed from Sydney to Parramatta to get all of the values in.
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Parramatta&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
## # A tibble: 1 x 2
##     lon   lat
##   <dbl> <dbl>
## 1  151. -33.8
map <- get_googlemap(center = c(151.0,-33.8), zoom = 10) #using the results from the code above, sub in your coordinates to set the centre of the map.
## Source : https://maps.googleapis.com/maps/api/staticmap?center=-33.8,151&zoom=10&size=640x640&scale=2&maptype=terrain&key=xxx-MtKgp2mleNTDL2te3Y-JNBSX2klxVY0
# ggmap has the same syntax as ggplot

ggmap(map) +
  geom_point(data = sc2, aes(x= lon, y = lat, size = Pop, color = IncomeValue)) #this line adds plot points. Set a variable to "size" and "colour" to map the values of that variable using size and colour of plots. Make sure this is put inside the aes() brackets.

ggmap(map) +
  geom_point(data = sc2, aes(x= lon, y = lat, size = Pop, color = IncomeValue))

ggmap(map) +
  geom_point(data = sc2, aes(x= lon, y = lat, shape = Population, color = Income), size = 3) +
  labs(caption = "(based on data collected by the Australian Bureau of Statistics in the 2016 census)", x = "Longitude", y = "Latitude")

str(sc2)
## tibble [36 × 9] (S3: tbl_df/tbl/data.frame)
##  $ SC         : chr [1:36] "Westfield Eastgardens" "Fashion Spree Orange Grove" "Neeta City" "Marrickville Metro" ...
##  $ Suburb     : chr [1:36] "Eastgardens" "Liverpool" "Fairfield" "Marrickville" ...
##  $ Person     : chr [1:36] "Tarryn" "Tarryn" "Tarryn" "Tarryn" ...
##  $ Population : chr [1:36] "LOW" "HIGH" "LOW" "HIGH" ...
##  $ Income     : chr [1:36] "HIGH" "LOW" "LOW" "HIGH" ...
##  $ Pop        : num [1:36] 860 27804 18081 26592 10468 ...
##  $ IncomeValue: num [1:36] 1833 1089 975 1814 2559 ...
##  $ lon        : num [1:36] 151 151 151 151 151 ...
##  $ lat        : num [1:36] -33.9 -33.9 -33.9 -33.9 -33.8 ...

mapview


library(ggplot2)




geocode("Parramatta") #check coordinates of where you want your centre of the map to be. I changed from Sydney to Parramatta to get all of the values in.

map <- get_googlemap(center = c(151.0,-33.8), zoom = 10) #using the results from the code above, sub in your coordinates to set the centre of the map.

# ggmap has the same syntax as ggplot

ggmap(map) +
  geom_point(data = sc2, aes(x= lon, y = lat, size = Pop, color = IncomeValue)) #this line adds plot points. Set a variable to "size" and "colour" to map the values of that variable using size and colour of plots. Make sure this is put inside the aes() brackets.
 

ggmap(map) +
  geom_point(data = sc2, aes(x= lon, y = lat, size = Pop, color = IncomeValue))

ggmap(map) +
  geom_point(data = sc2, aes(x= lon, y = lat, shape = Population, color = Income), size = 3) +
  labs(caption = "(based on data collected by the Australian Bureau of Statistics in the 2016 census)", x = "Longitude", y = "Latitude")






Habitat Variables

library(readr)
hab1 <- read_csv("withfactors.csv", col_types = cols(Population = col_factor(levels = c("HIGH", 
    "LOW")), Income = col_factor(levels = c("HIGH", 
    "LOW"))))
## Warning: Missing column names filled in: 'X16' [16], 'X17' [17], 'X18' [18],
## 'X19' [19]

PCA

library(REdaS)
## Loading required package: grid
cormatrix <- cor(hab1[1:36,4:15])
cormatrix
##                  carparks escalators   elevators      levels     toilets
## carparks       1.00000000  0.1104018 -0.03651974  0.27488040  0.36134990
## escalators     0.11040178  1.0000000  0.73954196  0.50882483  0.71726211
## elevators     -0.03651974  0.7395420  1.00000000  0.44497674  0.47146012
## levels         0.27488040  0.5088248  0.44497674  1.00000000  0.38411470
## toilets        0.36134990  0.7172621  0.47146012  0.38411470  1.00000000
## hours          0.12710142  0.1515900  0.09310969  0.01157244  0.08566331
## eateries       0.40250032  0.4862020  0.31493815  0.49847231  0.67115063
## trainstations  0.18535517  0.3602172  0.22379553  0.27026228  0.34410700
## busstops       0.12063463  0.1247951  0.15460860 -0.08987194 -0.02223400
## cinema         0.44525146  0.3371226  0.22596641  0.68586483  0.50876630
## greenspace    -0.05889480 -0.2018897 -0.07751053 -0.22812489 -0.14558848
## foodcourts     0.42506151  0.4548239  0.07409223  0.31467775  0.37917306
##                     hours   eateries trainstations    busstops      cinema
## carparks       0.12710142  0.4025003     0.1853552  0.12063463  0.44525146
## escalators     0.15159000  0.4862020     0.3602172  0.12479512  0.33712259
## elevators      0.09310969  0.3149381     0.2237955  0.15460860  0.22596641
## levels         0.01157244  0.4984723     0.2702623 -0.08987194  0.68586483
## toilets        0.08566331  0.6711506     0.3441070 -0.02223400  0.50876630
## hours          1.00000000  0.2923818     0.1756103  0.39344563  0.09786211
## eateries       0.29238185  1.0000000     0.2516909  0.11067248  0.69984069
## trainstations  0.17561033  0.2516909     1.0000000 -0.12707898  0.30303450
## busstops       0.39344563  0.1106725    -0.1270790  1.00000000 -0.14285714
## cinema         0.09786211  0.6998407     0.3030345 -0.14285714  1.00000000
## greenspace    -0.28375623 -0.3502501    -0.2115139 -0.17131295 -0.14151940
## foodcourts     0.10887017  0.4216912     0.3893209  0.06992765  0.34963827
##                greenspace foodcourts
## carparks      -0.05889480 0.42506151
## escalators    -0.20188969 0.45482390
## elevators     -0.07751053 0.07409223
## levels        -0.22812489 0.31467775
## toilets       -0.14558848 0.37917306
## hours         -0.28375623 0.10887017
## eateries      -0.35025006 0.42169121
## trainstations -0.21151389 0.38932091
## busstops      -0.17131295 0.06992765
## cinema        -0.14151940 0.34963827
## greenspace     1.00000000 0.06131806
## foodcourts     0.06131806 1.00000000
library(ggcorrplot)
ggcorrplot(cormatrix)

ggcorrplot_clustered <- ggcorrplot(cormatrix, hc.order = TRUE, type  ="lower")
ggcorrplot_clustered

Hierarchical clustering of correlation. Eelvators and escalators highly correlated, eateries and toilets, cinema and eateries and toilets and escaltors are all highly correlated.

bart_spher(cormatrix)
##  Bartlett's Test of Sphericity
## 
## Call: bart_spher(x = cormatrix)
## 
##      X2 = 325.279
##      df = 66
## p-value < 2.22e-16

Bartlett’s Test of Sphericity = testing assumptions to see if the data is correlated.

p-value = <2.222e<16

Make PCs

pca1 <- prcomp(cor(hab1[1:36,4:15]), scale = TRUE)
pca2 <- princomp(cor(hab1[1:36,4:15]), cor = TRUE)


summary(pca1)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.3153 1.4838 1.3814 0.96836 0.83908 0.59568 0.55923
## Proportion of Variance 0.4467 0.1835 0.1590 0.07814 0.05867 0.02957 0.02606
## Cumulative Proportion  0.4467 0.6302 0.7892 0.86734 0.92601 0.95558 0.98164
##                           PC8     PC9    PC10    PC11      PC12
## Standard deviation     0.3811 0.21125 0.15612 0.07828 6.375e-17
## Proportion of Variance 0.0121 0.00372 0.00203 0.00051 0.000e+00
## Cumulative Proportion  0.9937 0.99746 0.99949 1.00000 1.000e+00

PC1 explains 44.67% of the variance, PC2 explains 18.35% of the total variance, PC3 explains 15.90% and PC4 explains 7.81% of the total variance.These 4 PCs explain a cumulative proportion of 86.7% of the variance.

Screeplot

screeplot(pca1, type = "lines")

pca1$x
##                      PC1         PC2         PC3        PC4         PC5
## carparks       0.1359197  1.42762788 -2.12064115 -0.6480406  0.42151020
## escalators    -1.7758270 -0.98825800  1.47672596  0.4624595  0.84534876
## elevators     -0.3488096 -1.45542020  2.42897906 -0.3110332  0.01612957
## levels        -1.6302632  0.45123704  0.63126664 -0.8687229 -0.82116336
## toilets       -1.9955015  0.02737616  0.53731394 -0.1275671  0.62717405
## hours          1.9504040 -2.18421583 -1.54202171  0.4378214 -0.98924196
## eateries      -2.0217986 -0.31036094 -0.88292413 -0.7082315  0.02691283
## trainstations -0.4321372  0.40809526  0.08074933  2.2304057 -1.21981355
## busstops       3.3206029 -2.24103373 -0.84759577 -0.4980453  0.75193287
## cinema        -1.7808635  1.08700961 -0.49980569 -1.0038613 -0.92467426
## greenspace     5.1417437  2.39232169  1.65758265 -0.3054750 -0.10713926
## foodcourts    -0.5634696  1.38562107 -0.91962911  1.3402903  1.37302411
##                       PC6         PC7         PC8         PC9          PC10
## carparks      -0.24735751  1.07550188 -0.52082425 -0.06842893  0.0307474622
## escalators     0.02636193 -0.09866588 -0.25190649  0.16329447 -0.1002073999
## elevators      0.05055560  0.38441168 -0.29792874 -0.40997153  0.0003711169
## levels         1.13341732  0.01701822 -0.17031582  0.30597723  0.1837423642
## toilets       -1.15276440  0.07944638  0.16885678  0.31872145 -0.0163579956
## hours         -0.37532273 -0.70890312 -0.56631294  0.06364605 -0.0118396531
## eateries      -0.40709415 -0.55098661  0.46695602 -0.24158190  0.2843990468
## trainstations -0.04749158  0.67870695  0.37529164 -0.01684887  0.0341781632
## busstops       0.55132449  0.44108681  0.60034682  0.08590062 -0.0629123851
## cinema         0.10271361 -0.29838173  0.25561152 -0.10693748 -0.3663499310
## greenspace    -0.33440181 -0.35220951  0.03814829  0.02502054  0.0469579081
## foodcourts     0.70005923 -0.66702509 -0.09792283 -0.11879164 -0.0227286967
##                       PC11          PC12
## carparks       0.036258646  3.191891e-16
## escalators     0.178261906 -5.689893e-16
## elevators     -0.080050390 -1.665335e-16
## levels        -0.036670109  5.342948e-16
## toilets       -0.120354080  1.838807e-16
## hours         -0.018750742  3.330669e-16
## eateries       0.068548932 -6.383782e-16
## trainstations  0.017061656 -1.040834e-17
## busstops      -0.008841529 -5.273559e-16
## cinema         0.006878692  3.903128e-17
## greenspace     0.032259949  1.387779e-16
## foodcourts    -0.074602930  3.330669e-16
habpcscores <- cbind(hab1[1:36,1:15], pca1$x)
## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded
habpcscores
##                         Shopping Centre Population Income carparks escalators
## 1                 Westfield Eastgardens        LOW   HIGH        1          9
## 2            Fashion Spree Orange Grove       HIGH    LOW        3          0
## 3                            Neeta City        LOW    LOW        1         10
## 4                    Marrickville Metro       HIGH   HIGH        1          8
## 5                    Kellyville Village        LOW   HIGH        1          1
## 6                   Blacktown Westpoint       HIGH    LOW        5          9
## 7                     Penrith Westfield        LOW    LOW        6         19
## 8        Royal Randwick Shopping Centre       HIGH   HIGH        1          0
## 9                Rouse Hill Town Centre        LOW   HIGH        4          5
## 10                 Hurstville Westfield       HIGH    LOW       10         14
## 11               Revesby Village Centre        LOW    LOW        1          6
## 12          Bridgepoint Shopping Centre       HIGH   HIGH        1          2
## 13                     Birkenhead Point        LOW   HIGH        3         16
## 14            Westfield Shopping Centre       HIGH    LOW        1         21
## 15                    Minto Marketplace        LOW    LOW        2          2
## 16                        Castle Towers       HIGH   HIGH        3          1
## 17             Westfield Bondi Junction        LOW   HIGH        3         49
## 18                     Rhodes Waterside       HIGH    LOW        1          3
## 19                             Top Ryde        LOW    LOW        1          6
## 20                 Westfield Parramatta       HIGH   HIGH        2          4
## 21                         DFO Homebush        LOW   HIGH        2          8
## 22                    Bankstown Central       HIGH    LOW        2         24
## 23                    Westfield Burwood        LOW    LOW        5         21
## 24                        Ashfield Mall       HIGH   HIGH        1          6
## 25                     Macquarie Centre        LOW    LOW        5         17
## 26                       Warringah Mall        LOW   HIGH        6          7
## 27                     Westield Miranda       HIGH    LOW        3         28
## 28                        Dee Why Grand       HIGH   HIGH        1          4
## 29 Stockland Piccadilly Shopping Centre        LOW   HIGH        1          2
## 30                        South Village       HIGH    LOW        1         12
## 31               Westfield Mount Druitt        LOW    LOW        7          4
## 32                  Westfield Chatswood       HIGH   HIGH        1         64
## 33                   Kings Cross Centre        LOW   HIGH        1          0
## 34                       Auburn Central       HIGH    LOW        1          6
## 35             Broadway Shopping Centre        LOW    LOW        2          2
## 36                    Strathfield Plaza       HIGH   HIGH        1          1
##    elevators levels toilets hours eateries trainstations busstops cinema
## 1          6      3      21  12.0       52             1        1      1
## 2          0      1       6  11.0        6             0        1      0
## 3          3      3       5  12.0       19             1        1      0
## 4          2      2       2  12.0       20             0        1      0
## 5          2      1       1  12.0       12             0        1      0
## 6          7      9       6  12.0       32             1        1      1
## 7          4      5       8  12.0       45             1        1      1
## 8          2      2       6  12.0       19             0        1      0
## 9          2      4      10  12.0       31             1        1      1
## 10         3      6      14  11.5       48             1        1      1
## 11         1      3       4  12.0        3             1        1      0
## 12         2      2       2  12.0        6             0        1      0
## 13         5      4       6   9.5       18             0        1      0
## 14         4      4       5  12.0       65             1        1      1
## 15         3      1       4   8.5       16             1        1      0
## 16         1      4       7  12.0       61             0        1      1
## 17         4      7      24  11.5       47             1        1      1
## 18         1      7       3  10.0       38             0        1      1
## 19         6      7       6  12.0       63             0        1      1
## 20         3      6       3   9.5       26             1        1      1
## 21         3      2      12   8.0       15             0        1      0
## 22         5      3      17  12.0       46             1        1      1
## 23         5      5      23  11.5       41             1        1      1
## 24         6      4       7  12.0       18             1        1      0
## 25         5      5      24  11.5      104             1        1      1
## 26         6      3      13  12.0       44             0        1      1
## 27         8      5      22  11.5       74             1        1      1
## 28         2      2       6  11.0        8             0        1      0
## 29         4      2       3  11.0       10             1        1      0
## 30        10      3       4  11.5       11             1        1      0
## 31         3      3       9  12.0       47             1        1      1
## 32        29      8      24  11.5       57             1        1      1
## 33         2      3       1  12.0       16             1        1      1
## 34         1      3       3  12.0       14             1        1      0
## 35         6      5       5   9.0        7             0        1      1
## 36         0      5      10   8.5       17             1        0      1
##    greenspace foodcourts        PC1         PC2         PC3        PC4
## 1           0          1  0.1359197  1.42762788 -2.12064115 -0.6480406
## 2           0          1 -1.7758270 -0.98825800  1.47672596  0.4624595
## 3           0          1 -0.3488096 -1.45542020  2.42897906 -0.3110332
## 4           0          1 -1.6302632  0.45123704  0.63126664 -0.8687229
## 5           0          0 -1.9955015  0.02737616  0.53731394 -0.1275671
## 6           0          2  1.9504040 -2.18421583 -1.54202171  0.4378214
## 7           0          2 -2.0217986 -0.31036094 -0.88292413 -0.7082315
## 8           0          0 -0.4321372  0.40809526  0.08074933  2.2304057
## 9           1          2  3.3206029 -2.24103373 -0.84759577 -0.4980453
## 10          0          2 -1.7808635  1.08700961 -0.49980569 -1.0038613
## 11          0          1  5.1417437  2.39232169  1.65758265 -0.3054750
## 12          3          1 -0.5634696  1.38562107 -0.91962911  1.3402903
## 13          1          2  0.1359197  1.42762788 -2.12064115 -0.6480406
## 14          0          3 -1.7758270 -0.98825800  1.47672596  0.4624595
## 15          1          2 -0.3488096 -1.45542020  2.42897906 -0.3110332
## 16          0          1 -1.6302632  0.45123704  0.63126664 -0.8687229
## 17          0          3 -1.9955015  0.02737616  0.53731394 -0.1275671
## 18          0          1  1.9504040 -2.18421583 -1.54202171  0.4378214
## 19          0          1 -2.0217986 -0.31036094 -0.88292413 -0.7082315
## 20          0          1 -0.4321372  0.40809526  0.08074933  2.2304057
## 21          1          1  3.3206029 -2.24103373 -0.84759577 -0.4980453
## 22          0          2 -1.7808635  1.08700961 -0.49980569 -1.0038613
## 23          1          2  5.1417437  2.39232169  1.65758265 -0.3054750
## 24          0          1 -0.5634696  1.38562107 -0.91962911  1.3402903
## 25          0          2  0.1359197  1.42762788 -2.12064115 -0.6480406
## 26          1          2 -1.7758270 -0.98825800  1.47672596  0.4624595
## 27          0          1 -0.3488096 -1.45542020  2.42897906 -0.3110332
## 28          0          1 -1.6302632  0.45123704  0.63126664 -0.8687229
## 29          0          0 -1.9955015  0.02737616  0.53731394 -0.1275671
## 30          1          2  1.9504040 -2.18421583 -1.54202171  0.4378214
## 31          0          1 -2.0217986 -0.31036094 -0.88292413 -0.7082315
## 32          0          1 -0.4321372  0.40809526  0.08074933  2.2304057
## 33          1          1  3.3206029 -2.24103373 -0.84759577 -0.4980453
## 34          0          1 -1.7808635  1.08700961 -0.49980569 -1.0038613
## 35          1          0  5.1417437  2.39232169  1.65758265 -0.3054750
## 36          1          1 -0.5634696  1.38562107 -0.91962911  1.3402903
##            PC5         PC6         PC7         PC8         PC9          PC10
## 1   0.42151020 -0.24735751  1.07550188 -0.52082425 -0.06842893  0.0307474622
## 2   0.84534876  0.02636193 -0.09866588 -0.25190649  0.16329447 -0.1002073999
## 3   0.01612957  0.05055560  0.38441168 -0.29792874 -0.40997153  0.0003711169
## 4  -0.82116336  1.13341732  0.01701822 -0.17031582  0.30597723  0.1837423642
## 5   0.62717405 -1.15276440  0.07944638  0.16885678  0.31872145 -0.0163579956
## 6  -0.98924196 -0.37532273 -0.70890312 -0.56631294  0.06364605 -0.0118396531
## 7   0.02691283 -0.40709415 -0.55098661  0.46695602 -0.24158190  0.2843990468
## 8  -1.21981355 -0.04749158  0.67870695  0.37529164 -0.01684887  0.0341781632
## 9   0.75193287  0.55132449  0.44108681  0.60034682  0.08590062 -0.0629123851
## 10 -0.92467426  0.10271361 -0.29838173  0.25561152 -0.10693748 -0.3663499310
## 11 -0.10713926 -0.33440181 -0.35220951  0.03814829  0.02502054  0.0469579081
## 12  1.37302411  0.70005923 -0.66702509 -0.09792283 -0.11879164 -0.0227286967
## 13  0.42151020 -0.24735751  1.07550188 -0.52082425 -0.06842893  0.0307474622
## 14  0.84534876  0.02636193 -0.09866588 -0.25190649  0.16329447 -0.1002073999
## 15  0.01612957  0.05055560  0.38441168 -0.29792874 -0.40997153  0.0003711169
## 16 -0.82116336  1.13341732  0.01701822 -0.17031582  0.30597723  0.1837423642
## 17  0.62717405 -1.15276440  0.07944638  0.16885678  0.31872145 -0.0163579956
## 18 -0.98924196 -0.37532273 -0.70890312 -0.56631294  0.06364605 -0.0118396531
## 19  0.02691283 -0.40709415 -0.55098661  0.46695602 -0.24158190  0.2843990468
## 20 -1.21981355 -0.04749158  0.67870695  0.37529164 -0.01684887  0.0341781632
## 21  0.75193287  0.55132449  0.44108681  0.60034682  0.08590062 -0.0629123851
## 22 -0.92467426  0.10271361 -0.29838173  0.25561152 -0.10693748 -0.3663499310
## 23 -0.10713926 -0.33440181 -0.35220951  0.03814829  0.02502054  0.0469579081
## 24  1.37302411  0.70005923 -0.66702509 -0.09792283 -0.11879164 -0.0227286967
## 25  0.42151020 -0.24735751  1.07550188 -0.52082425 -0.06842893  0.0307474622
## 26  0.84534876  0.02636193 -0.09866588 -0.25190649  0.16329447 -0.1002073999
## 27  0.01612957  0.05055560  0.38441168 -0.29792874 -0.40997153  0.0003711169
## 28 -0.82116336  1.13341732  0.01701822 -0.17031582  0.30597723  0.1837423642
## 29  0.62717405 -1.15276440  0.07944638  0.16885678  0.31872145 -0.0163579956
## 30 -0.98924196 -0.37532273 -0.70890312 -0.56631294  0.06364605 -0.0118396531
## 31  0.02691283 -0.40709415 -0.55098661  0.46695602 -0.24158190  0.2843990468
## 32 -1.21981355 -0.04749158  0.67870695  0.37529164 -0.01684887  0.0341781632
## 33  0.75193287  0.55132449  0.44108681  0.60034682  0.08590062 -0.0629123851
## 34 -0.92467426  0.10271361 -0.29838173  0.25561152 -0.10693748 -0.3663499310
## 35 -0.10713926 -0.33440181 -0.35220951  0.03814829  0.02502054  0.0469579081
## 36  1.37302411  0.70005923 -0.66702509 -0.09792283 -0.11879164 -0.0227286967
##            PC11          PC12
## 1   0.036258646  3.191891e-16
## 2   0.178261906 -5.689893e-16
## 3  -0.080050390 -1.665335e-16
## 4  -0.036670109  5.342948e-16
## 5  -0.120354080  1.838807e-16
## 6  -0.018750742  3.330669e-16
## 7   0.068548932 -6.383782e-16
## 8   0.017061656 -1.040834e-17
## 9  -0.008841529 -5.273559e-16
## 10  0.006878692  3.903128e-17
## 11  0.032259949  1.387779e-16
## 12 -0.074602930  3.330669e-16
## 13  0.036258646  3.191891e-16
## 14  0.178261906 -5.689893e-16
## 15 -0.080050390 -1.665335e-16
## 16 -0.036670109  5.342948e-16
## 17 -0.120354080  1.838807e-16
## 18 -0.018750742  3.330669e-16
## 19  0.068548932 -6.383782e-16
## 20  0.017061656 -1.040834e-17
## 21 -0.008841529 -5.273559e-16
## 22  0.006878692  3.903128e-17
## 23  0.032259949  1.387779e-16
## 24 -0.074602930  3.330669e-16
## 25  0.036258646  3.191891e-16
## 26  0.178261906 -5.689893e-16
## 27 -0.080050390 -1.665335e-16
## 28 -0.036670109  5.342948e-16
## 29 -0.120354080  1.838807e-16
## 30 -0.018750742  3.330669e-16
## 31  0.068548932 -6.383782e-16
## 32  0.017061656 -1.040834e-17
## 33 -0.008841529 -5.273559e-16
## 34  0.006878692  3.903128e-17
## 35  0.032259949  1.387779e-16
## 36 -0.074602930  3.330669e-16
pc1aov <- aov(habpcscores$PC1 ~ busstops*greenspace*hours, data = habpcscores)
summary(pc1aov)
##                  Df Sum Sq Mean Sq F value Pr(>F)  
## busstops          1   0.33   0.327   0.072 0.7906  
## greenspace        1  24.18  24.179   5.309 0.0281 *
## hours             1   8.71   8.707   1.912 0.1766  
## greenspace:hours  1   2.51   2.509   0.551 0.4635  
## Residuals        31 141.17   4.554                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pc2aov <- aov(habpcscores$PC2 ~ carparks*greenspace*cinema, data = habpcscores)
summary(pc2aov)
##                            Df Sum Sq Mean Sq F value Pr(>F)
## carparks                    1   0.00  0.0008   0.000  0.986
## greenspace                  1   0.01  0.0118   0.005  0.946
## cinema                      1   0.08  0.0750   0.030  0.864
## carparks:greenspace         1   0.02  0.0210   0.008  0.928
## carparks:cinema             1   1.46  1.4556   0.579  0.453
## greenspace:cinema           1   0.25  0.2531   0.101  0.753
## carparks:greenspace:cinema  1   0.44  0.4435   0.176  0.678
## Residuals                  28  70.39  2.5140
pc3aov <- aov(habpcscores$PC3 ~ escalators*elevators, data = habpcscores)
summary(pc3aov)
##                      Df Sum Sq Mean Sq F value Pr(>F)
## escalators            1   0.19  0.1935   0.102  0.751
## elevators             1   1.12  1.1213   0.591  0.448
## escalators:elevators  1   0.98  0.9845   0.519  0.476
## Residuals            32  60.67  1.8960
pc4aov <- aov(habpcscores$PC4 ~ foodcourts*trainstations*hours, data = habpcscores)
summary(pc4aov)
##                                Df Sum Sq Mean Sq F value Pr(>F)  
## foodcourts                      1  1.168   1.168   1.414 0.2444  
## trainstations                   1  0.195   0.195   0.236 0.6307  
## hours                           1  0.469   0.469   0.568 0.4575  
## foodcourts:trainstations        1  0.127   0.127   0.154 0.6978  
## foodcourts:hours                1  0.003   0.003   0.003 0.9545  
## trainstations:hours             1  3.892   3.892   4.714 0.0385 *
## foodcourts:trainstations:hours  1  1.973   1.973   2.390 0.1333  
## Residuals                      28 23.118   0.826                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In PC1, greenspace is significant (p = 0.0281).

In PC4, we see an interaction between presence of train stations and opening hours on Thursdays with a p-value of 0.0385.

library(FactoMineR)
sc_pca <- PCA(hab1 [1:36,4:15])

summary(sc_pca)
## 
## Call:
## PCA(X = hab1[1:36, 4:15]) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5   Dim.6   Dim.7
## Variance               4.351   1.583   1.422   1.087   0.958   0.650   0.618
## % of var.             36.262  13.190  11.852   9.057   7.986   5.416   5.151
## Cumulative % of var.  36.262  49.452  61.304  70.360  78.347  83.763  88.913
##                        Dim.8   Dim.9  Dim.10  Dim.11  Dim.12
## Variance               0.520   0.366   0.214   0.150   0.080
## % of var.              4.337   3.049   1.784   1.249   0.668
## Cumulative % of var.  93.251  96.299  98.084  99.332 100.000
## 
## Individuals (the 10 first)
##                   Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
## 1             |  2.573 |  1.242  0.985  0.233 |  0.601  0.633  0.054 | -0.274
## 2             |  2.991 | -2.424  3.751  0.657 |  0.661  0.767  0.049 |  0.581
## 3             |  2.092 | -1.045  0.698  0.250 |  1.084  2.062  0.269 | -0.335
## 4             |  2.668 | -1.992  2.533  0.558 |  1.303  2.981  0.239 | -0.096
## 5             |  3.548 | -2.932  5.487  0.683 |  1.506  3.980  0.180 | -0.242
## 6             |  3.358 |  2.011  2.582  0.359 | -0.148  0.038  0.002 |  0.796
## 7             |  2.562 |  1.863  2.217  0.529 |  0.184  0.059  0.005 |  1.356
## 8             |  3.154 | -2.400  3.678  0.579 |  1.378  3.332  0.191 | -0.287
## 9             |  2.138 |  0.589  0.222  0.076 | -0.443  0.344  0.043 |  1.343
## 10            |  4.112 |  2.589  4.277  0.396 | -0.405  0.289  0.010 |  2.419
##                  ctr   cos2  
## 1              0.147  0.011 |
## 2              0.659  0.038 |
## 3              0.219  0.026 |
## 4              0.018  0.001 |
## 5              0.115  0.005 |
## 6              1.237  0.056 |
## 7              3.593  0.280 |
## 8              0.161  0.008 |
## 9              3.524  0.395 |
## 10            11.429  0.346 |
## 
## Variables (the 10 first)
##                  Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## carparks      |  0.484  5.385  0.234 | -0.101  0.641  0.010 |  0.641 28.847
## escalators    |  0.785 14.173  0.617 |  0.103  0.674  0.011 | -0.472 15.633
## elevators     |  0.569  7.453  0.324 |  0.152  1.453  0.023 | -0.693 33.766
## levels        |  0.718 11.848  0.516 | -0.229  3.325  0.053 | -0.104  0.763
## toilets       |  0.800 14.698  0.640 | -0.089  0.505  0.008 | -0.138  1.348
## hours         |  0.269  1.662  0.072 |  0.711 31.918  0.505 |  0.259  4.710
## eateries      |  0.820 15.468  0.673 |  0.081  0.414  0.007 |  0.199  2.787
## trainstations |  0.517  6.149  0.268 | -0.080  0.406  0.006 |  0.025  0.043
## busstops      |  0.084  0.163  0.007 |  0.788 39.209  0.621 |  0.099  0.684
## cinema        |  0.747 12.818  0.558 | -0.302  5.771  0.091 |  0.246  4.263
##                 cos2  
## carparks       0.410 |
## escalators     0.222 |
## elevators      0.480 |
## levels         0.011 |
## toilets        0.019 |
## hours          0.067 |
## eateries       0.040 |
## trainstations  0.001 |
## busstops       0.010 |
## cinema         0.061 |
sc_pca$eig
##         eigenvalue percentage of variance cumulative percentage of variance
## comp 1  4.35145114             36.2620928                          36.26209
## comp 2  1.58280047             13.1900039                          49.45210
## comp 3  1.42218382             11.8515319                          61.30363
## comp 4  1.08682219              9.0568516                          70.36048
## comp 5  0.95832816              7.9860680                          78.34655
## comp 6  0.64994119              5.4161765                          83.76272
## comp 7  0.61808282              5.1506902                          88.91341
## comp 8  0.52045800              4.3371500                          93.25056
## comp 9  0.36583013              3.0485844                          96.29915
## comp 10 0.21413594              1.7844662                          98.08362
## comp 11 0.14983211              1.2486009                          99.33222
## comp 12 0.08013402              0.6677835                         100.00000

Varimax

fa1 <- factanal(hab1[1:36,4:15], 4, rotation = "varimax")
fa1
## 
## Call:
## factanal(x = hab1[1:36, 4:15], factors = 4, rotation = "varimax")
## 
## Uniquenesses:
##      carparks    escalators     elevators        levels       toilets 
##         0.651         0.005         0.357         0.394         0.346 
##         hours      eateries trainstations      busstops        cinema 
##         0.791         0.151         0.785         0.795         0.070 
##    greenspace    foodcourts 
##         0.760         0.005 
## 
## Loadings:
##               Factor1 Factor2 Factor3 Factor4
## carparks               0.438   0.357   0.151 
## escalators     0.924   0.114   0.297   0.199 
## elevators      0.788                   0.104 
## levels         0.454   0.606   0.144  -0.113 
## toilets        0.587   0.394   0.208   0.334 
## hours                                  0.444 
## eateries       0.254   0.672   0.236   0.527 
## trainstations  0.259   0.212   0.321         
## busstops              -0.172           0.409 
## cinema         0.206   0.930   0.150         
## greenspace    -0.153  -0.163   0.141  -0.412 
## foodcourts     0.146   0.189   0.967         
## 
##                Factor1 Factor2 Factor3 Factor4
## SS loadings      2.254   2.193   1.430   1.014
## Proportion Var   0.188   0.183   0.119   0.084
## Cumulative Var   0.188   0.371   0.490   0.574
## 
## Test of the hypothesis that 4 factors are sufficient.
## The chi square statistic is 26.01 on 24 degrees of freedom.
## The p-value is 0.352

nMDS

library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-6
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(ggplot2)

library(readr)
full <- read_csv("full.csv", col_types = cols(Population = col_factor(levels = c("LOW", 
    "HIGH")), Income = col_factor(levels = c("LOW", 
    "HIGH")), PopulationValue = col_number(), 
    IncomeValue = col_number()))


library(readr)
assemblage <- read_csv("assemblage.csv")
## Warning: Missing column names filled in: 'X38' [38], 'X39' [39]
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_double(),
##   Species = col_character(),
##   X38 = col_logical(),
##   X39 = col_logical()
## )
## ℹ Use `spec()` for the full column specifications.
library(readr)
species <- read_csv("species.csv")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_double(),
##   ShoppingCentre = col_character(),
##   Suburb = col_character(),
##   Population = col_character(),
##   Income = col_character(),
##   PopulationValue = col_number(),
##   IncomeValue = col_number()
## )
## ℹ Use `spec()` for the full column specifications.
spp1 <- (na.omit( species
        [2:32,4:20])) #square brackets used to select only numeric data

library(readr)
spc <- read_csv("spp.csv")
## Warning: Missing column names filled in: 'X19' [19], 'X20' [20]
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_double(),
##   `Shopping Centre` = col_character(),
##   Suburb = col_character(),
##   Population = col_character(),
##   Income = col_character(),
##   X19 = col_logical(),
##   X20 = col_logical()
## )
## ℹ Use `spec()` for the full column specifications.
spp <- spc[1:36,5:18]

library(readr)
factors <- read_csv("factors.csv", col_types = cols(Population = col_factor(levels = c("LOW", 
    "HIGH")), Income = col_factor(levels = c("LOW", 
    "HIGH"))))
## Warning: Missing column names filled in: 'X1' [1], 'X7' [7], 'X8' [8], 'X9' [9],
## 'X10' [10], 'X11' [11], 'X12' [12], 'X13' [13], 'X14' [14], 'X15' [15],
## 'X16' [16], 'X17' [17], 'X18' [18], 'X19' [19], 'X20' [20], 'X21' [21],
## 'X22' [22], 'X23' [23], 'X24' [24], 'X25' [25], 'X26' [26], 'X27' [27],
## 'X28' [28], 'X29' [29], 'X30' [30], 'X31' [31], 'X32' [32], 'X33' [33],
## 'X34' [34], 'X35' [35], 'X36' [36], 'X37' [37], 'X38' [38], 'X39' [39]
str(spp)
## tibble [36 × 14] (S3: tbl_df/tbl/data.frame)
##  $ Department: num [1:36] 4 0 0 2 0 5 3 0 2 2 ...
##  $ Grocery   : num [1:36] 5 0 4 6 3 10 3 3 5 7 ...
##  $ Accessory : num [1:36] 22 14 6 5 1 23 30 3 22 21 ...
##  $ Clothing  : num [1:36] 36 31 15 7 1 37 55 9 34 30 ...
##  $ Pharmacy  : num [1:36] 4 0 1 1 1 2 3 1 2 3 ...
##  $ Pet       : num [1:36] 1 0 0 0 0 0 2 0 0 0 ...
##  $ Toy       : num [1:36] 2 0 0 1 0 1 1 0 3 1 ...
##  $ Specialty : num [1:36] 7 0 1 2 4 6 25 2 2 1 ...
##  $ Thrift    : num [1:36] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Services  : num [1:36] 63 2 26 32 18 63 62 22 36 44 ...
##  $ Homeware  : num [1:36] 6 5 0 1 0 6 4 1 8 3 ...
##  $ Technology: num [1:36] 8 1 4 6 1 10 10 1 9 9 ...
##  $ Cosmetic  : num [1:36] 2 0 0 0 0 5 4 0 5 8 ...
##  $ Fresh Food: num [1:36] 52 6 19 20 15 53 70 19 6 43 ...
TransSpp <- spp^(1/4)
spp.dis <- vegdist(TransSpp, method ="bray")
spp.mds0 <- isoMDS(spp.dis)
## initial  value 13.795483 
## iter   5 value 10.457427
## iter  10 value 10.286952
## iter  15 value 9.983493
## final  value 9.873226 
## converged

Plotting nMDS

plot <- ordiplot(spp.mds0, type = "t")
## species scores not available

pchs <- c(0:5)
Pop_Factor <- factors$Population
nMDS_spp<- metaMDS(TransSpp, distance="bray", k=2)
## Run 0 stress 0.0969794 
## Run 1 stress 0.126783 
## Run 2 stress 0.1072847 
## Run 3 stress 0.1071067 
## Run 4 stress 0.105539 
## Run 5 stress 0.1250816 
## Run 6 stress 0.1342471 
## Run 7 stress 0.1001295 
## Run 8 stress 0.1066917 
## Run 9 stress 0.1062038 
## Run 10 stress 0.1062039 
## Run 11 stress 0.1066749 
## Run 12 stress 0.1357998 
## Run 13 stress 0.1199346 
## Run 14 stress 0.1198358 
## Run 15 stress 0.1274054 
## Run 16 stress 0.1171052 
## Run 17 stress 0.1166742 
## Run 18 stress 0.1072848 
## Run 19 stress 0.1062039 
## Run 20 stress 0.1267099 
## *** No convergence -- monoMDS stopping criteria:
##     20: stress ratio > sratmax
plot_spp <- ordiplot(nMDS_spp, display = "sites", type = "n")
points(nMDS_spp, col = "black", pch = pchs[Pop_Factor])
legend("topright", bty = "n", legend = levels(Pop_Factor), pch = pchs)

nMDS_spp
## 
## Call:
## metaMDS(comm = TransSpp, distance = "bray", k = 2) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     TransSpp 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.0969794 
## Stress type 1, weak ties
## No convergent solutions - best solution after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'TransSpp'
plot_spp1 <- ordiplot(nMDS_spp, display = "species", type = "n")
points(nMDS_spp, col = "black", pch = pchs[Pop_Factor])
legend("topright", bty = "n", legend = levels(Pop_Factor), pch = pchs)

nMDS_spp
## 
## Call:
## metaMDS(comm = TransSpp, distance = "bray", k = 2) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     TransSpp 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.0969794 
## Stress type 1, weak ties
## No convergent solutions - best solution after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'TransSpp'

Income as a Factor

plotinc <- ordiplot(spp.mds0, type = "t")
## species scores not available

pchs <- c(0:5)
Income_Factor <- factors$Income
nMDS_spp<- metaMDS(TransSpp, distance="bray", k=2)
## Run 0 stress 0.0969794 
## Run 1 stress 0.1266484 
## Run 2 stress 0.1266776 
## Run 3 stress 0.1236921 
## Run 4 stress 0.1270373 
## Run 5 stress 0.09910931 
## Run 6 stress 0.1236998 
## Run 7 stress 0.134555 
## Run 8 stress 0.1066722 
## Run 9 stress 0.1198341 
## Run 10 stress 0.1266631 
## Run 11 stress 0.09911034 
## Run 12 stress 0.1284552 
## Run 13 stress 0.1051754 
## Run 14 stress 0.1260333 
## Run 15 stress 0.09978155 
## Run 16 stress 0.1076054 
## Run 17 stress 0.1075795 
## Run 18 stress 0.1062512 
## Run 19 stress 0.1055391 
## Run 20 stress 0.09697946 
## ... Procrustes: rmse 3.860572e-05  max resid 0.0001006442 
## ... Similar to previous best
## *** Solution reached
plot_sppincome <- ordiplot(nMDS_spp, display = "species", type = "n")
points(nMDS_spp, col = "black", pch = pchs[Income_Factor])
legend("topright", bty = "n", legend = levels(Income_Factor), pch = pchs)

PERMANOVA

adonis(TransSpp ~ Income*Population, data = factors, permutations = 999)
## 
## Call:
## adonis(formula = TransSpp ~ Income * Population, data = factors,      permutations = 999) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                   Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## Income             1   0.02536 0.025357 0.50286 0.01492  0.697
## Population         1   0.02312 0.023118 0.45845 0.01360  0.744
## Income:Population  1   0.03767 0.037667 0.74697 0.02216  0.516
## Residuals         32   1.61366 0.050427         0.94932       
## Total             35   1.69980                  1.00000

No significant difference between groups.

Species Diversity

library(readr)
simp <- read_csv("simp.csv", col_types = cols(Population = col_factor(levels = c("LOW", 
    "HIGH")), Income = col_factor(levels = c("LOW", 
    "HIGH"))))


```r
shapiro.test(simp$SimpsonsD)
## 
##  Shapiro-Wilk normality test
## 
## data:  simp$SimpsonsD
## W = 0.90166, p-value = 0.003789
shapiro.test(simp$Eveness)
## 
##  Shapiro-Wilk normality test
## 
## data:  simp$Eveness
## W = 0.96847, p-value = 0.3852
simp.aov <- aov(SimpsonsD ~ Income*Population, data = simp)
summary(simp.aov)
##                   Df  Sum Sq  Mean Sq F value Pr(>F)
## Income             1 0.00208 0.002079   0.805  0.376
## Population         1 0.00052 0.000520   0.201  0.657
## Income:Population  1 0.00741 0.007413   2.869  0.100
## Residuals         32 0.08270 0.002584

No significant interaction or main effect on diversity.