# Calculate distance (in miles) based on latitude and longitude of locations dist = function(lat1, long1, lat2, long2){ # R = 3963; R = 6378.7; Dlong = long2 - long1; a = sin(lat1/57.3)*sin(lat2/57.3) + cos(lat1/57.3)*cos(lat2/57.3)*cos(Dlong/57.3); d = R*acos(a); d; } # read data dl = read.table("diversity.txt", header=T) intro = read.table("introduction.txt", header=T) # Sites located within radius are considered as neighbors radius = 160 n1 = dim(dl)[1] n2 = dim(intro)[1] size = rep(0, n1) nevent = rep(0, n1) sizemat = matrix(0, nrow=n1, ncol=4) neventmat = matrix(0, nrow=n1, ncol=4) neventnewmat = matrix(0, nrow=n1, ncol=4) dismat = matrix(0, nrow=n1, ncol=n2) # Calculate total introduction sizes and introduction effort for neighbors for (i in 1:n1){ for (j in 1:n2){ dis = dist(dl$Latitude[i], dl$Longitude[i], intro$Latitude[j], intro$Longitude[j]) dismat[i,j] = dis if (dis <= radius){ size[i] = size[i] + intro$Size[j] id = intro$YearID[j] sizemat[i,id] = sizemat[i,id] + intro$Size[j] } } } radius1= 390 # Calculate total introduction sizes and introduction effort for neighbors for (i in 1:n1){ for (j in 1:n2){ # dis = dist(dl$Latitude[i], dl$Longitude[i], intro$Latitude[j], intro$Longitude[j]) # dismat[i,j] = dis dis = dismat[i,j] if (dis <= radius1){ nevent[i] = nevent[i] + 1 id = intro$YearID[j] neventmat[i,id] = neventmat[i,id] + 1 } } } dnear = rep(0, n1) for (i in 1:n1){ dnear[i] = min(dismat[i,]) } dl$size = size dl$dnear = dnear dl$nevent = nevent #SIZE# library(betareg) size10 = log10(size+1) cf = coef(lm(qlogis(diversity)~size10,data=dl)) fm2 = betareg(diversity~size10, data=dl, link="logit", start=list(mean=cf, precision=1)) summary(fm2) #SIZEMAT# sizemat=sizemat+1 sizemat[,1]=log10(sizemat[,1]) sizemat[,2]=log10(sizemat[,2]) sizemat[,3]=log10(sizemat[,3]) sizemat[,4]=log10(sizemat[,4]) fm21 = betareg(diversity ~ sizemat[,1]+sizemat[,2]+sizemat[,3]+sizemat[,4],data=dl, link="logit") summary(fm21) #dnear# dnear20=dnear+1 fm31 = betareg(diversity ~ dnear20, data=dl, link="logit") summary(fm31) #nevent# neventnew = log10(nevent+1) # You can choose whatever transformation you want cf = coef(lm(qlogis(diversity) ~ neventnew, data = dl)) fm4 = betareg(diversity ~ neventnew, data=dl, link="logit", start=list(mean=cf, precision=1)) summary(fm4) neventnewmat[,1] = log10(neventmat[,1]+1) neventnewmat[,2] = log10(neventmat[,2]+1) neventnewmat[,3] = log10(neventmat[,3]+1) neventnewmat[,4] = log10(neventmat[,4]+1) fm41 = betareg(diversity ~ neventnewmat[,1]+neventnewmat[,2]+neventnewmat[,3]+neventnewmat[,4],data=dl, link="logit") summary(fm41) #################################### # R code for pearson's correlation # #################################### corPerm1 <- function(x,y,nperm) # # En - This function computes a two-tailed permutation test of a Pearson # correlation coefficient between two data vectors. The test statistic is r. # # Fr - Cette fonction realise, par permutation, un test bilateral du coefficient # de correlation de Pearson entre deux vecteurs. La statistique du test est r. # # Parameters of the function: # x, y: the two data vectors # nperm = number of permutations # # Example: test the correlation between two vectors of random numbers N(0,1) # x <- rnorm(50,0,1) # y <- rnorm(50,0,1) # cor.out = corPerm1(x, y, 999) # cor.out # Compare the result to: cor.test(x,y) # # Pedagogic objectives -- # - illustrate the structure of a function in the R language # - show how to program a "for" loop (command: "for(i in 1:nperm)") # - show how to write out the results using "cat" # - show how to program a permutation test # # Pierre Legendre, October 2005 { x <- as.matrix(x) y <- as.matrix(y) n <- nrow(x) r.ref <- cor(x,y) nGT <- 1 for(i in 1:nperm) { y.perm <- sample(y,n) r.perm <- cor(x,y.perm) if( abs(r.perm) >= abs(r.ref) ) nGT <- nGT+1 } P <- nGT/(nperm+1) cat('\nPearson correlation (two-tailed test)','\n') cat('r =',r.ref,'\n') cat('Prob(',nperm,'permutations) =',P,'\n','\n') return(list( Correlation=r.ref, No.perm=nperm, P.perm=P )) } corPerm1(dl$diversity, size10, nperm=1000) corPerm1(dl$diversity, neventnew, nperm=1000) corPerm1(dl$diversity, dnear20, nperm=1000) ################################################################################## # specify neighbors for each sampling site based on its longitude and latitude. # The weight matrix is based on inverse distance ################################################################################## library(spdep) coords = as.matrix(cbind(dl$Longitude, dl$Latitude)) coord.nb = dnearneigh(coords, 0, 4100, longlat=TRUE) dsts = nbdists(coord.nb, coords, longlat=T) idw = lapply(dsts, function(x) 1/x) dist.list = nb2listw(coord.nb, glist=idw) # Moran'I of the spatial autocorrelation in original data # test based on randomization moran.test(dl$diversity, dist.list, alternative="two.sided",zero.policy=T) # test based on permutation moran.mc(dl$diversity, dist.list, zero.policy=T, nsim=1000) ################################################################################# # correlation between h and distance from each sampling site to Meadowlands, NJ # ################################################################################# disvec = rep(0,n1) for (i in 1:n1){ disvec[i] = dist(dl$Latitude[i], dl$Longitude[i], 40.818138, -74.07390) } corPerm1(dl$diversity, disvec, nperm=1000)