This is my first time trying to plot an RDA in R. My datas for the RDA are two matrix; one with environmental conditions such as pH, temperature, etc. and the other with different species (relative frequencies of every species at each station). These data were taken at different stations. My response variable is the species and my explanatory variable is the environmental conditions.
> chalut_large_mat #Species found at the stations
Actiniaria sp. Chionoecetes opilio Lycodes sp. Meganyctiphanes norvegica Pandalus borealis
EM14 0.04347826 0.1304348 0.04347826 0.04347826 0.6521739
EM16 0.00000000 0.0629275 0.00000000 0.00000000 0.8440492
BIC3 0.00000000 0.0000000 0.00000000 0.00000000 0.8431373
EM17 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
EM12 0.00000000 0.4076087 0.00000000 0.00000000 0.0000000
Pseudopleuronectes americanus Aspidophoroides monopterygius Crevette sp. Clymenella torquata Eteone sp.
EM14 0.08695652 0.00000000 0.00000000 0.0000000 0.00000000
EM16 0.00000000 0.03146375 0.00000000 0.0000000 0.00000000
BIC3 0.00000000 0.00000000 0.05882353 0.0000000 0.00000000
EM17 0.00000000 0.00000000 0.00000000 0.2142857 0.07142857
EM12 0.04347826 0.00000000 0.00000000 0.0000000 0.00000000
Mya arenaria Nuculana sp. Gammarus sp. Echinoidea Zoarcidae Sclerocrangon boreasÿ Pandalus montagui
EM14 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.0000000
EM16 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.0000000
BIC3 0.0000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.0000000
EM17 0.5714286 0.07142857 0.07142857 0.00000000 0.00000000 0.00000000 0.0000000
EM12 0.0000000 0.00000000 0.00000000 0.06521739 0.05978261 0.03804348 0.1521739
Sebastes sp.
EM14 0.0000000
EM16 0.0000000
BIC3 0.0000000
EM17 0.0000000
EM12 0.1684783
> chimie_large_mat #environmental conditions measured at every station
prof groupe lat long temp sal pH MPS chla phaeo oxygen_abs DO
BIC3 293 1 48.61018 -68.95022 5.9218 34.4832 7.885 11.29020 0.08353328 0.010738708 0.7734 53.06538
EM12 115 3 48.81703 -68.75740 3.3890 33.4904 7.986 15.95486 0.06623374 0.012311763 1.3739 135.49446
EM14 327 3 48.70100 -68.65102 6.2420 34.7703 7.919 15.72575 0.04919374 0.013765077 0.7829 49.63710
EM16 72 4 48.58262 -68.53828 7.4710 27.6365 8.197 9.96400 0.25115478 0.009154732 2.0803 284.46747
EM17 22 4 48.54998 -68.56083 6.6148 28.1667 8.207 11.71675 0.59121391 0.004247401 2.1196 270.89014
DIC AT phosphate NO2.NO3 pCO2
BIC3 2541.000 2275.107 2.481480 25.163658 1023.3
EM12 2380.742 2222.264 1.809933 18.547943 758.3
EM14 2494.135 2255.000 2.523323 24.916785 924.3
EM16 2297.125 2139.078 1.113962 7.995123 461.3
EM17 2221.665 2067.097 1.129026 10.197057 432.8
> ####ACR####
> acr<-rda(chalut_large_mat,chimie_large_mat,scale=FALSE)
> summary(acr)
Call:
rda(X = chalut_large_mat, Y = chimie_large_mat, scale = FALSE)
Partitioning of variance:
Inertia Proportion
Total 0.3105 1
Constrained 0.3105 1
Unconstrained 0.0000 0
Eigenvalues, and their contribution to the variance
Importance of components:
RDA1 RDA2 RDA3 RDA4
Eigenvalue 0.2307 0.07607 0.002925 0.0008345
Proportion Explained 0.7429 0.24499 0.009419 0.0026876
Cumulative Proportion 0.7429 0.98789 0.997312 1.0000000
Accumulated constrained eigenvalues
Importance of components:
RDA1 RDA2 RDA3 RDA4
Eigenvalue 0.2307 0.07607 0.002925 0.0008345
Proportion Explained 0.7429 0.24499 0.009419 0.0026876
Cumulative Proportion 0.7429 0.98789 0.997312 1.0000000
Scaling 2 for species and site scores
* Species are scaled proportional to eigenvalues
* Sites are unscaled: weighted dispersion equal on all dimensions
* General scaling constant of scores: 1.055675
Species scores
RDA1 RDA2 RDA3 RDA4
Actiniaria sp. -0.0094665 0.003555 0.0350172 -0.0053410
Chionoecetes opilio 0.1196407 0.296348 0.0175344 0.0253551
Lycodes sp. -0.0094665 0.003555 0.0350172 -0.0053410
Meganyctiphanes norvegica -0.0094665 0.003555 0.0350172 -0.0053410
Pandalus borealis -0.8158654 -0.104382 -0.0029265 0.0111654
Pseudopleuronectes americanus -0.0004895 0.038002 0.0621369 -0.0110708
Aspidophoroides monopterygius -0.0124340 -0.001963 -0.0067394 0.0225114
Crevette sp. -0.0237428 -0.007329 -0.0242368 -0.0357637
Clymenella torquata 0.1269299 -0.129706 0.0005282 0.0052062
Eteone sp. 0.0423100 -0.043235 0.0001761 0.0017354
Mya arenaria 0.3384797 -0.345882 0.0014086 0.0138833
Nuculana sp. 0.0423100 -0.043235 0.0001761 0.0017354
Gammarus sp. 0.0423100 -0.043235 0.0001761 0.0017354
Echinoidea 0.0276653 0.046338 -0.0118461 -0.0005831
Zoarcidae 0.0253599 0.042476 -0.0108590 -0.0005345
Sclerocrangon boreasÿ 0.0161381 0.027030 -0.0069102 -0.0003402
Pandalus montagui 0.0645524 0.108121 -0.0276410 -0.0013607
Sebastes sp. 0.0714687 0.119705 -0.0306025 -0.0015064
Site scores (weighted sums of species scores)
RDA1 RDA2 RDA3 RDA4
EM14 -0.2426 0.09112 0.897572 -0.136903
EM16 -0.4404 -0.06954 -0.238709 0.797357
BIC3 -0.4498 -0.13884 -0.459181 -0.677566
EM17 0.6601 -0.67457 0.002747 0.027076
EM12 0.4728 0.79183 -0.202429 -0.009965
Site constraints (linear combinations of constraining variables)
RDA1 RDA2 RDA3 RDA4
EM14 -0.2426 0.09112 0.897572 -0.136903
EM16 -0.4404 -0.06954 -0.238709 0.797357
BIC3 -0.4498 -0.13884 -0.459181 -0.677566
EM17 0.6601 -0.67457 0.002747 0.027076
EM12 0.4728 0.79183 -0.202429 -0.009965
Biplot scores for constraining variables
RDA1 RDA2 RDA3 RDA4
prof -0.7333 -0.20114 0.2821 -0.5850
groupe 0.6258 -0.02513 -0.7714 0.1125
lat -0.7837 -0.24602 -0.3465 0.4530
long 0.6651 -0.06890 -0.7343 -0.1173
I tried biplot but I get this error message and I know that this is not the way I want to represent my data.
>biplot(acr, scaling = 1)
Error in biplot.rda(acr, scaling = 1, main = "RDA – scaling 1") :
'biplot.rda' not suitable for models with constraints
For the plot, I want both my stations and my species as points and all my environmental conditions as vectors. I know I can do this with ggplot but I have no idea how. I tried but I only get my stations as points and my species as vector. I don't know how to put the environmental conditions as vector. This is the code I did to create my plot but it is obviously wrong. Thank you for your help!
acpl_summ<-summary(acr,scaling=1)
acpl_sites<-acpl_summ$sites[,1:2]
head(acpl_sites)
sites<-as_tibble(acpl_sites)%>%
mutate(station=rownames(acpl_sites))
sites
##Vector##
acpl_sp<-acpl_summ$species[,1:2]
acpl_sp
especes<-as_tibble(acpl_sp)%>%
mutate(espece=rownames(acpl_sp))
especes
##Plot##
acp_graph_sites<-ggplot(sites)+
geom_text(aes(x=RDA1,y=RDA2,label=station))+
xlab('RDA 1 (xx%)')+
ylab('RDA 2(xx%)')+
geom_hline(yintercept=0,linetype='dotted')+
geom_vline(xintercept=0,linetype='dotted')+
theme_bw()+
theme(panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
legend.position='top',
legend.text=element_text(size=7))+
scale_colour_discrete(name=NULL)
acp_graph_sites
exp_sp<-4 #expansion factor
acp_graph_sites_sp<-acp_graph_sites+
geom_text(data=especes,aes(x=RDA1/exp_sp,
y=RDA2/exp_sp,
label=espece))+
geom_segment(data=especes,
aes(x=0,y=0,xend=RDA1/exp_sp,yend=RDA2/exp_sp))
acp_graph_sites_sp
