Lovell — Oct 14, 2014, 8:34 AM
rm(list=ls())
setwd("~/Desktop/Adaptomics_Project/niche_modeling")
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.0-10
library(fields)
Loading required package: spam
Warning: package 'spam' was built under R version 3.1.1
Loading required package: grid
Spam version 1.0-1 (2014-09-09) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
Attaching package: 'spam'
The following objects are masked from 'package:base':
backsolve, forwardsolve
Loading required package: maps
Warning: package 'maps' was built under R version 3.1.1
library(varSelRF)
Loading required package: randomForest
Warning: package 'randomForest' was built under R version 3.1.1
randomForest 4.6-10
Type rfNews() to see new features/changes/bug fixes.
#read in the data
acc.data<-read.csv("mr2014_accessioninfo2.csv",header=T)
acc.data<-acc.data[,c("lat","lon","elev",paste("bio",c(1:19), sep=""),"apollo","taxon","ploidy"),]
for (i in unique(acc.data$taxon)){
data.in<-acc.data[acc.data$taxon==i,]
tab.in<-table(data.in$apollo, data.in$ploidy)
if(dim(tab.in)[2]>0){
print(i)
print(tab.in)
}
}
[1] "Boechera stricta"
2 3
0 106 0
1 14 2
[1] "Boechera pinetorum"
2 3
0 1 3
1 0 13
[1] "Boechera retrofracta"
2 3
0 30 0
1 27 35
[1] "Boechera divaricarpa"
2 3 4
0 0 9 2
1 3 26 3
[1] "Boechera divaricarpa x Boechera retrofracta"
3
1 1
[1] "Boechera species"
2 3
0 4 1
1 3 1
[1] "Boechera collinsii"
2 3
0 3 2
1 0 1
[1] "Boechera pendulocarpa"
2 3
0 4 1
1 1 0
[1] "Boechera stricta x Boechera retrofracta"
2
0 1
1 9
[1] "Boechera fernaldiana ssp. fernaldiana"
2 3
0 6 0
1 1 1
[1] "Boechera breweri"
2
0 2
1 0
[1] "Boechera cobrensis"
2
0 2
1 0
[1] "Boechera crandallii"
2
0 20
1 4
[1] "Boechera lemmonii"
2 3
0 3 0
1 0 1
[1] "Boechera lignifera"
2 3
0 1 0
1 3 3
[1] "Boechera sparsiflora"
2 3
0 0 2
1 1 0
[1] "Boechera microphylla"
2
0 0
1 4
[1] "Boechera pallidifolia"
2 3
0 10 0
1 12 3
[1] "Boechera pendulina"
2
0 4
1 1
[1] "Boechera perennans"
2
0 3
1 1
[1] "Boechera fendleri"
2 3
0 1 1
1 3 0
[1] "Boechera puberula"
2 3
0 3 0
1 0 1
[1] "Boechera williamsii"
2
0 5
1 13
[1] "Boechera laevigata"
2
0 6
1 0
[1] "Boechera missouriensis"
2
0 1
[1] "Boechera lyallii"
2
0 2
1 2
[1] "Boechera canadensis"
2
0 7
1 0
[1] "Boechera platysperma"
2
0 2
1 0
[1] "Boechera rectissima"
2
0 1
1 1
[1] "Boechera schistacea"
2
0 1
[1] "Boechera shockleyi"
2
0 2
[1] "Boechera shortii"
2
0 1
[1] "Boechera suffrutescens"
2
0 0
1 0
[1] "Polyctenium fremontii var. confertum"
3
0 1
[1] "Boechera inyoensis"
3
0 0
1 1
[1] "Boechera lincolnensis"
3
0 1
1 1
[1] "Boechera formosa"
2
0 5
1 1
[1] "Boechera xylopoda"
2
0 1
[1] "Boechera glaucovalvula"
3
0 1
[1] "Boechera gracilipes"
2
0 1
1 1
[1] "Boechera patens"
2
0 1
[1] "B. oxylobula"
2
0 1
[1] "Boechera falcifructa"
3
1 1
[1] "Boechera glareosa"
2
1 1
[1] "Boechera lasiocarpa"
2
0 9
[1] "Boechera johnstonii"
2
0 1
[1] "Boechera nevadensis"
2
0 1
[1] "Boechera oxylobula"
2
0 10
[1] "Boechera polyantha"
2
0 5
1 2
[1] "Boechera polyantha x Boechera retrofracta"
3
1 2
[1] "Boechera stricta x Boechera spatifolia"
2 3
1 10 5
[1] "Boechera spatifolia"
2 3
0 8 0
1 8 3
#first for apollo-
#make sure data is complete for apollo
apollo.data<-acc.data[,-which(names(acc.data)=="ploidy")]
apollo.data<-apollo.data[complete.cases(apollo.data[,-which(names(apollo.data) %in% c("taxon")),]),]
write.csv(apollo.data,"mr2014_apolloniche.csv",row.names=F)
#data subset
apollo.space<-apollo.data[,c("lat","lon")]
apollo.env<-apollo.data[,c("elev",paste("bio",c(1:19), sep=""))]
apollo.gen<-apollo.data[,c("apollo")]
#random forest variable selection- all species
full.forest<-randomForest(x=apollo.env, y=as.factor(apollo.gen))
varImpPlot(full.forest, main="Full Sample Random Forest Variable Importance")
apollo.rfsel<-varSelRF(xdata=apollo.env, Class=as.factor(apollo.gen),vars.drop.frac = 0.2)
apollo.rfsel
Backwards elimination on random forest; ntree = 5000 ; mtryFactor = 1
Selected variables:
[1] "bio1" "bio10" "bio15" "elev"
Number of selected variables: 4
sel.vars<-apollo.rfsel$selected.vars
print(sel.vars)
[1] "bio1" "bio10" "bio15" "elev"
best.env<-apollo.env[,sel.vars]
#cca multivariate test w/o geo
cca1<-cca(best.env~apollo.gen)
aov.out.apollo<-anova(cca1,by="terms",permu=10000)
print(aov.out.apollo)
Permutation test for cca under reduced model
Terms added sequentially (first to last)
Model: cca(formula = best.env ~ apollo.gen)
Df Chisq F N.Perm Pr(>F)
apollo.gen 1 0.00 0.87 9999 0.075 .
Residual 1593 0.06
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#cca multivariate test w/ geo
geo.dist<-rdist.earth(apollo.space[,c("lat","lon")], miles=F)
cca2<-cca(best.env~as.factor(apollo.gen),z=geo.dist)
aov.out.apollo.partial<-anova(cca2,by="terms",permu=10000)
print(aov.out.apollo.partial)
Permutation test for cca under reduced model
Terms added sequentially (first to last)
Model: cca(formula = best.env ~ as.factor(apollo.gen), z = geo.dist)
Df Chisq F N.Perm Pr(>F)
as.factor(apollo.gen) 1 0.00 0.87 9999 0.081 .
Residual 1593 0.06
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#get the species to loop through
#criteria here is we need 10 observations of each apollo 1 and apollo 2
species.to.test<-vector()
for (i in unique(apollo.data$taxon)){
data.in<-apollo.data[apollo.data$taxon==i,]
tab.out<-as.numeric(table(data.in$apollo))
sex<-tab.out[1]
apo<-tab.out[2]
if(is.na(sex) | is.na(apo)){next
}else{
if(sex>=5 & apo>=5) {species.to.test<-c(species.to.test,i)
}else{next}
}
}
print(species.to.test)
[1] "Boechera stricta" "Boechera pinetorum"
[3] "Boechera retrofracta" "Boechera divaricarpa"
[5] "Boechera species" "Boechera collinsii"
[7] "Boechera pendulocarpa" "Boechera crandallii"
[9] "Boechera lemmonii" "Boechera pauciflora"
[11] "Boechera sparsiflora" "Boechera microphylla"
[13] "Boechera pallidifolia" "Boechera perennans"
[15] "Boechera fendleri" "Boechera puberula"
[17] "Boechera williamsii" "Boechera lyallii"
[19] "Boechera spatifolia"
species.to.test<-species.to.test[-which(species.to.test=="Boechera microphylla")]
for (i in species.to.test){
apollo.data.in<-apollo.data[apollo.data$taxon==i,]
apollo.space.in<-apollo.data.in[,c("lat","lon")]
apollo.env.in<-apollo.data.in[,c("elev",paste("bio",c(1:19), sep=""))]
apollo.gen.in<-apollo.data.in[,c("apollo")]
apollo.rfsel.in<-varSelRF(xdata=apollo.env.in, Class=as.factor(apollo.gen.in),vars.drop.frac = 0.2)
full.forest.in<-randomForest(x=apollo.env, y=as.factor(apollo.gen))
#varImpPlot(full.forest.in, main=paste(i,"Random Forest Variable Importance"))
sel.vars.in<-apollo.rfsel.in$selected.vars
best.env.in<-apollo.env.in[,sel.vars.in]
cat("for", i , "the selected environmental variables are:",sel.vars.in, "\n")
cca1<-cca(best.env.in~apollo.gen.in)
aov.out.apollo<-anova(cca1,by="terms",permu=10000)
f.raw<-data.frame(aov.out.apollo)$F[1]
pval.raw<-data.frame(aov.out.apollo)$Pr..F.[1]
geo.dist<-rdist.earth(apollo.space.in[,c("lat","lon")], miles=F)
cca2<-cca(best.env.in~as.factor(apollo.gen.in),z=geo.dist)
aov.out.apollo.partial<-anova(cca2,by="terms",permu=10000)
f.part<-data.frame(aov.out.apollo.partial)$F[1]
pval.part<-data.frame(aov.out.apollo.partial)$Pr..F.[1]
cat("n.apo0 = ",length(apollo.gen.in[apollo.gen.in==0]), "... n.apo1 = ",length(apollo.gen.in[apollo.gen.in==1]),"\n")
cat("f=",f.raw,"p.w/o geo = ", pval.raw,"f.part=",f.part, "p.partial geo = ", pval.part)
cat("\n*********\n")
}
for Boechera stricta the selected environmental variables are: bio16 bio18
n.apo0 = 182 ... n.apo1 = 32
f= 0.1092 p.w/o geo = 0.8697 f.part= 0.1092 p.partial geo = 0.8721
*********
for Boechera pinetorum the selected environmental variables are: bio14 bio17 bio3
n.apo0 = 20 ... n.apo1 = 57
f= 0.1339 p.w/o geo = 0.7956 f.part= 0.1339 p.partial geo = 0.7947
*********
for Boechera retrofracta the selected environmental variables are: bio3 elev
n.apo0 = 97 ... n.apo1 = 158
f= 19.8 p.w/o geo = 1e-04 f.part= 19.8 p.partial geo = 1e-04
*********
for Boechera divaricarpa the selected environmental variables are: bio14 bio17 bio18 bio4
n.apo0 = 46 ... n.apo1 = 158
f= 3.389 p.w/o geo = 0.0667 f.part= 3.389 p.partial geo = 0.0705
*********
for Boechera species the selected environmental variables are: bio10 bio5
n.apo0 = 8 ... n.apo1 = 19
f= 2.728 p.w/o geo = 0.1143 f.part= 2.728 p.partial geo = 0.1126
*********
for Boechera collinsii the selected environmental variables are: bio3 elev
n.apo0 = 8 ... n.apo1 = 9
f= 0.3736 p.w/o geo = 0.298 f.part= 0.3736 p.partial geo = 0.2942
*********
for Boechera pendulocarpa the selected environmental variables are: bio16 bio19 elev
n.apo0 = 32 ... n.apo1 = 8
f= 2.313 p.w/o geo = 0.1096 f.part= 2.313 p.partial geo = 0.1137
*********
for Boechera crandallii the selected environmental variables are: bio16 bio3
n.apo0 = 26 ... n.apo1 = 7
f= 6.173 p.w/o geo = 0.0197 f.part= 6.173 p.partial geo = 0.0208
*********
for Boechera lemmonii the selected environmental variables are: bio5 bio7 elev
n.apo0 = 23 ... n.apo1 = 7
f= 2.221 p.w/o geo = 0.0136 f.part= 2.221 p.partial geo = 0.0124
*********
for Boechera pauciflora the selected environmental variables are: bio10 bio17 bio4 bio9
n.apo0 = 6 ... n.apo1 = 18
f= 1.312 p.w/o geo = 0.2786 f.part= 1.312 p.partial geo = 0.2788
*********
for Boechera sparsiflora the selected environmental variables are: bio12 bio18 bio19 bio5 bio8
n.apo0 = 25 ... n.apo1 = 20
f= 0.7525 p.w/o geo = 0.4623 f.part= 0.7525 p.partial geo = 0.465
*********
for Boechera pallidifolia the selected environmental variables are: bio5 bio8
n.apo0 = 16 ... n.apo1 = 19
f= 2.08 p.w/o geo = 0.1382 f.part= 2.08 p.partial geo = 0.1358
*********
for Boechera perennans the selected environmental variables are: bio18 bio3
n.apo0 = 20 ... n.apo1 = 11
f= 5.875 p.w/o geo = 0.0396 f.part= 5.875 p.partial geo = 0.0341
*********
for Boechera fendleri the selected environmental variables are: bio3 bio9
n.apo0 = 8 ... n.apo1 = 8
f= 2.56 p.w/o geo = 0.1146 f.part= 2.56 p.partial geo = 0.1124
*********
for Boechera puberula the selected environmental variables are: bio5 bio9
n.apo0 = 30 ... n.apo1 = 17
f= 1.823 p.w/o geo = 0.1303 f.part= 1.823 p.partial geo = 0.1303
*********
for Boechera williamsii the selected environmental variables are: bio2 elev
n.apo0 = 8 ... n.apo1 = 15
f= 13.9 p.w/o geo = 0.001 f.part= 13.9 p.partial geo = 9e-04
*********
for Boechera lyallii the selected environmental variables are: bio18 bio8
n.apo0 = 14 ... n.apo1 = 5
f= NA p.w/o geo = NA f.part= NA p.partial geo = NA
*********
for Boechera spatifolia the selected environmental variables are: bio15 bio17
n.apo0 = 8 ... n.apo1 = 11
f= 8.191 p.w/o geo = 0.0117 f.part= 8.191 p.partial geo = 0.0124
*********
#make sure data is complete for ploidy
ploidy.data<-acc.data[,-which(names(acc.data)=="apollo")]
ploidy.data<-ploidy.data[complete.cases(ploidy.data[,-which(names(ploidy.data) %in% c("taxon")),]),]
write.csv(ploidy.data,"mr2014_ploidyniche.csv",row.names=F)
ploidy.space<-ploidy.data[,c("lat","lon")]
ploidy.env<-ploidy.data[,c("elev",paste("bio",c(1:19), sep=""))]
ploidy.gen<-ploidy.data[,c("ploidy")]
full.forest<-randomForest(x=ploidy.env, y=as.factor(ploidy.gen))
varImpPlot(full.forest, main="Full Sample Random Forest Variable Importance")
ploidy.rfsel<-varSelRF(xdata=ploidy.env, Class=as.factor(ploidy.gen),vars.drop.frac = 0.2)
ploidy.rfsel
Backwards elimination on random forest; ntree = 5000 ; mtryFactor = 1
Selected variables:
[1] "bio18" "bio4"
Number of selected variables: 2
sel.vars<-ploidy.rfsel$selected.vars
print(sel.vars)
[1] "bio18" "bio4"
best.env<-ploidy.env[,sel.vars]
cca1<-cca(best.env~ploidy.gen)
aov.out.ploidy<-anova(cca1,by="terms",permu=10000)
print(aov.out.ploidy)
Permutation test for cca under reduced model
Terms added sequentially (first to last)
Model: cca(formula = best.env ~ ploidy.gen)
Df Chisq F N.Perm Pr(>F)
ploidy.gen 1 0.00 0.39 9999 0.54
Residual 540 0.02
geo.dist<-rdist.earth(ploidy.space[,c("lat","lon")], miles=F)
cca2<-cca(best.env~as.factor(ploidy.gen),z=geo.dist)
aov.out.ploidy.partial<-anova(cca2,by="terms",permu=10000)
print(aov.out.ploidy.partial)
Permutation test for cca under reduced model
Terms added sequentially (first to last)
Model: cca(formula = best.env ~ as.factor(ploidy.gen), z = geo.dist)
Df Chisq F N.Perm Pr(>F)
as.factor(ploidy.gen) 2 0.00 2.10 9999 0.13
Residual 539 0.02
#get the species to loop through
#criteria here is we need 10 observations of each ploidy 1 and ploidy 2
species.to.test<-vector()
for (i in unique(ploidy.data$taxon)){
data.in<-ploidy.data[ploidy.data$taxon==i,]
tab.out<-as.numeric(table(data.in$ploidy))
sex<-tab.out[1]
apo<-tab.out[2]
if(is.na(sex) | is.na(apo)){next
}else{
if(sex>=3 & apo>=3) {species.to.test<-c(species.to.test,i)
}else{next}
}
}
print(species.to.test)
[1] "Boechera retrofracta"
[2] "Boechera collinsii"
[3] "Boechera divaricarpa"
[4] "Boechera lignifera"
[5] "Boechera pallidifolia"
[6] "Boechera stricta x Boechera spatifolia"
[7] "Boechera spatifolia"
species.to.test<-species.to.test[-which(species.to.test=="Boechera collinsii")]
#just retrofracta
for (i in species.to.test){
ploidy.data.in<-ploidy.data[ploidy.data$taxon==i,]
ploidy.space.in<-ploidy.data.in[,c("lat","lon")]
ploidy.env.in<-ploidy.data.in[,c("elev",paste("bio",c(1:19), sep=""))]
ploidy.gen.in<-ploidy.data.in[,c("ploidy")]
ploidy.rfsel.in<-varSelRF(xdata=ploidy.env.in, Class=as.factor(ploidy.gen.in),vars.drop.frac = 0.2)
full.forest.in<-randomForest(x=ploidy.env, y=as.factor(ploidy.gen))
#varImpPlot(full.forest.in, main=paste(i,"Random Forest Variable Importance"))
sel.vars.in<-ploidy.rfsel.in$selected.vars
best.env.in<-ploidy.env.in[,sel.vars.in]
cat("for", i , "the selected environmental variables are:",sel.vars.in, "\n")
cca1<-cca(best.env.in~ploidy.gen.in)
aov.out.ploidy<-anova(cca1,by="terms",permu=10000)
f.raw<-data.frame(aov.out.ploidy)$F[1]
pval.raw<-data.frame(aov.out.ploidy)$Pr..F.[1]
geo.dist<-rdist.earth(ploidy.space.in[,c("lat","lon")], miles=F)
cca2<-cca(best.env.in~as.factor(ploidy.gen.in),z=geo.dist)
aov.out.ploidy.partial<-anova(cca2,by="terms",permu=10000)
pval.part<-data.frame(aov.out.ploidy.partial)$Pr..F.[1]
f.part<-data.frame(aov.out.ploidy.partial)$F[1]
cat("n.2x = ",length(ploidy.gen.in[ploidy.gen.in==2]), "... n.3x = ",length(ploidy.gen.in[ploidy.gen.in==3]),"\n")
cat("f=",f.raw,"p.w/o geo = ", pval.raw,"f.part=",f.part, "p.partial geo = ", pval.part)
cat("\n*********\n")
}
for Boechera retrofracta the selected environmental variables are: bio3 elev
n.2x = 57 ... n.3x = 35
f= 13.58 p.w/o geo = 1e-04 f.part= 13.58 p.partial geo = 1e-04
*********
for Boechera divaricarpa the selected environmental variables are: bio19 bio4 bio9
n.2x = 3 ... n.3x = 35
f= 2.48 p.w/o geo = 0.1082 f.part= 7.212 p.partial geo = 5e-04
*********
for Boechera lignifera the selected environmental variables are: bio13 bio7
n.2x = 4 ... n.3x = 3
f= 25.1 p.w/o geo = 0.0311 f.part= 25.1 p.partial geo = 0.0301
*********
for Boechera pallidifolia the selected environmental variables are: bio15 bio17
n.2x = 22 ... n.3x = 3
f= 10.99 p.w/o geo = 1e-04 f.part= 10.99 p.partial geo = 5e-04
*********
for Boechera stricta x Boechera spatifolia the selected environmental variables are: bio1 bio14 bio17 bio2 bio5 bio9
n.2x = 10 ... n.3x = 5
f= 50.83 p.w/o geo = 0.0013 f.part= 50.83 p.partial geo = 9e-04
*********
for Boechera spatifolia the selected environmental variables are: bio16 bio4
n.2x = 16 ... n.3x = 3
f= 10.7 p.w/o geo = 0.0013 f.part= 10.7 p.partial geo = 0.0015
*********
#full analysis for just retrofracta
data.retro<-acc.data[acc.data$taxon=="Boechera retrofracta",]
data.retro<-data.retro[complete.cases(data.retro),]
write.csv(data.retro,"mr2014_retroniche.csv",row.names=F)
all.space.in<-data.retro[,c("lat","lon")]
all.env.in<-data.retro[,c("elev",paste("bio",c(1:19), sep=""))]
retro.ploidy<-data.retro[,c("ploidy")]
all.rfsel.in<-varSelRF(xdata=all.env.in, Class=as.factor(retro.ploidy),vars.drop.frac = 0.2)
full.forest.in<-randomForest(x=all.env.in, y=as.factor(retro.ploidy))
sel.vars.ploidy<-all.rfsel.in$selected.vars
retro.apollo<-data.retro[,c("apollo")]
all.rfsel.in<-varSelRF(xdata=all.env.in, Class=as.factor(retro.apollo),vars.drop.frac = 0.2)
full.forest.in<-randomForest(x=all.env.in, y=as.factor(retro.apollo))
sel.vars.apollo<-all.rfsel.in$selected.vars
retro.vars<-unique(c(sel.vars.apollo,sel.vars.ploidy))
retro.vars.data <- cbind(data.retro[,sel.vars.apollo])
ret.unique<-as.character(paste(data.retro$ploidy,data.retro$apollo, sep="."))
data.retro$ploidy.apollo<-ret.unique
retro.toplot<-data.retro[,c(retro.vars,"ploidy.apollo")]
library(GGally)
Warning: package 'GGally' was built under R version 3.1.1
ggpairs(retro.toplot, colour="ploidy.apollo",columns=c("elev","bio2","bio3"))
#manova
fit <- with(data.retro, manova(cbind(elev,bio2,bio3) ~ retro.ploidy*retro.apollo))
summary(fit, test="Pillai")
Df Pillai approx F num Df den Df Pr(>F)
retro.ploidy 1 0.402 19.5 3 87 9.2e-10 ***
retro.apollo 1 0.288 11.8 3 87 1.6e-06 ***
Residuals 89
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#just for diploids:
dip.retro<-data.retro[data.retro$ploidy==2,]
dip.retro<-dip.retro[complete.cases(dip.retro),]
all.space.in<-dip.retro[,c("lat","lon")]
all.env.in<-dip.retro[,c("elev",paste("bio",c(1:19), sep=""))]
retro.apollo<-dip.retro[,c("apollo")]
all.rfsel.in<-varSelRF(xdata=all.env.in, Class=as.factor(retro.apollo),vars.drop.frac = 0.2)
full.forest.in<-randomForest(x=all.env.in, y=as.factor(retro.apollo))
sel.vars.apollo<-all.rfsel.in$selected.vars
retro.toplot <- cbind(dip.retro[,sel.vars.apollo])
ret.apollo<-as.character(dip.retro$apollo)
retro.toplot$apollo<-ret.apollo
library(GGally)
ggpairs(retro.toplot, colour=as.character("apollo"),columns=sel.vars.apollo)
#manova
fit2 <- with(retro.toplot, manova(cbind(elev,bio2,bio3) ~ apollo))
summary(fit2, test="Pillai")
Df Pillai approx F num Df den Df Pr(>F)
apollo 1 0.321 8.36 3 53 0.00012 ***
Residuals 55
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit2 <- with(retro.toplot, manova(cbind(elev,bio2) ~ apollo))
summary(fit2, test="Pillai")
Df Pillai approx F num Df den Df Pr(>F)
apollo 1 0.314 12.4 2 54 3.8e-05 ***
Residuals 55
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1