In atredennick/community_synchrony: Calculates syunchrony and stability metrics from time series data

Test of INLA regression fitting functions

library(communitySynchrony)
library(INLA)
library(reshape2)
library(ggplot2)

Format data

sppList=sort(c("PSSP","HECO","POSE","ARTR"))
alpha.effect=c(0.014,0.017,0.026,0.018) # for spp in alphabetical order

# Set up list for storage
all_dataframes <- list()
all_crowding <- list()
for(spp in 1:length(sppList)){
  doSpp=sppList[spp]
  growDfile=paste("../../data/Idaho/",doSpp,"/growDnoNA.csv",sep="")
  growD=read.csv(growDfile)
  D <- growD  #subset(growD,allEdge==0)
  D$logarea.t0 <- log(D$area.t0)
  D$logarea.t1 <- log(D$area.t1)
  D$quad <- as.character(D$quad)

  # remove outliers (large plants that obviously do not turn into tiny plants) for ARTR only
  if(doSpp=="ARTR"){
    tmp=which(D$quad=="Q23" & D$year==45 & D$trackID==67)
    tmp=c(tmp,which(D$quad=="Q12" & D$year==55 & D$trackID==25))
    tmp=c(tmp,which(D$quad=="Q26" & D$year==45 & D$trackID==73))
    D=D[-tmp,]
  }else{
    D=D
  }

  # calculate crowding 
  for(i in 1:length(sppList)){
    distDfile=paste("../../data/Idaho/",sppList[i],"/",sppList[i],"_genet_xy.csv",sep="")
    if(i==1){
      distD=read.csv(distDfile)
      distD$nbSpp=sppList[i]  
    }else{
      tmp=read.csv(distDfile)
      tmp$nbSpp=sppList[i] 
      distD=rbind(distD,tmp)
    }
  }

  distD=distD[,c("quad","year","trackID","area","nbSpp","x","y")]
  W=matrix(NA,dim(D)[1],length(sppList))
  for(i in 1:dim(D)[1]){
    tmpD=subset(distD,year==D$year[i] & quad==D$quad[i])
    focal=which(tmpD$trackID==D$trackID[i] & tmpD$nbSpp==doSpp)
    xx=tmpD$x[focal] ; yy=tmpD$y[focal]
    tmpD$distance=sqrt((xx-tmpD$x)^2+(yy-tmpD$y)^2)
    tmpD=subset(tmpD,distance>0)
    if(dim(tmpD)[1]>0){
      for(k in 1:length(sppList)){
        sppI=which(tmpD$nbSpp==sppList[k])
        if(length(sppI)>0){
          W[i,k]=sum(exp(-1*alpha.effect[k]*tmpD$distance[sppI]^2)*tmpD$area[sppI])         
        }else{
          W[i,k]=0
        }
      }
    }else{
      W[i,]=0
    }   
  }

  all_crowding[[spp]] <- W
  all_dataframes[[spp]] <- D
}
names(all_crowding) <- sppList
names(all_dataframes) <- sppList

Call regression function for growth

growth_params <- list()
for(spp in sppList){
  tmp <- get_growth_params(dataframe = all_dataframes[[spp]],
                           crowd_mat = all_crowding[[spp]])
  growth_params[[spp]] <- tmp
}

Plot yearly intercepts and slopes

intercepts <- matrix(ncol=length(sppList), nrow=nrow(growth_params[[1]]))
slopes <- matrix(ncol=length(sppList), nrow=nrow(growth_params[[1]]))
i <- 1
for(spp in sppList){
  intercepts[,i] <- growth_params[[spp]]$Intercept.yr
  slopes[,i] <- growth_params[[spp]]$logarea.t0.yr
  i <- i+1
}
colnames(intercepts) <- sppList
colnames(slopes) <- sppList
intm <- melt(intercepts)
slopem <- melt(slopes)
intm$variable <- "intercept"
slopem$variable <- "slope"
combo <- rbind(intm, slopem)
ggplot(subset(combo, Var2!="ARTR"), aes(x=Var1, y=value, color=Var2))+
  geom_line()+
  geom_point()+
  facet_grid(variable~., scales="free_y")

atredennick/community_synchrony documentation built on May 10, 2019, 2:10 p.m.