データ・Rコード

発表論文で使用したデータやRコードを以下のサイトで公開しています。​

1. Chen et al. 2023: Nature (figshare)
2. Evans et al. 2022: Ecology Letters (figshare)
3. Hisano et al. 2021: Global Change Biology (figshare)
4. Hisano & Chen 2020: Global Ecology and Biogeography (figshare)
5. Hisano et al. 2019: Ecology Letters (figshare)

---

データ編集用Rコード・メモ

連番を追加する（data.table, dplyr）
org<-org[,":="(Census=row_number(Year_Mon),
                        Censuses=length(Year_Mon)),by=PLOT_ID]
​
文字列を取り出す
​soil[,":="(Month=substr(date,5,6))]

left関数、right関数
left<-function(x,chr){return(substr(x,1,chr))}​
​right<-function(x,chr){return(substr(x,nchar(x)-chr+1,nchar(x)))}

相互作用の簡易確認
library(sjPlot);library(sjmisc)
plot_model(m1_MEM, type = "pred", terms = c("Year", "CMIave"))

相互作用の簡易図示
interact_plot(mb, pred = Y, modx = logCVba, interval = TRUE)

mb<-lmer(AGBGI_MEMcorr~Year+sqrtSA1+log(SpCVba1+1)+
           Year:log(SpCVba1+1)+
           (1|PLOT_ID),can)

p<-interact_plot(mb,pred = Year,modx = SpCVba1,data=can,
              interval = TRUE,vary.lty = F,
              colors = c("salmon1","darkturquoise","mediumpurple1"),
              line.thickness = 1.5,
              y.label = expression(paste("Growth (Mg ", " ", ha^{-1}, " ", yr^{-1},")")),
              x.label = expression(paste("Calendar year")),
              legend.main = expression(paste(CV[BA])),
              modx.values = c(quantile(arid$SpCVba1,0.05),mean(arid$SpCVba1),quantile(arid$SpCVba1,0.95)),
              modx.labels = c("5th percentile","Mean","95th percentile"),
              main.title = expression(paste("Temporal trends dependent on ",CV[BA]))
              )+
  geom_line(aes(x=Year,y=AGBGI_MEMcorr),linetype = "solid",size=1.5)+
  theme_classic()
p​

カテゴリカル変数からダミー変数作成
library(fastDummies)
dat<-can[,.(Eco,Prov)]
dat1<-dummy_cols(dat)
can<-cbind(can,dat1)

変数の順序を変える
to_graph$red<-factor(to_graph$red, levels=c("Red","Non-red"))

data.frame（symbol）全体から複数の文字列（wik$SPEC）を検索する
which(Reduce(`|`, lapply(wik$SPEC, `==`, symbol)), arr.ind = TRUE)
もしくは、
symbol %>%
  mutate_all(list(~ . %in% wik$SPEC)) %>%
  as.matrix %>% 
  which(., arr.ind = TRUE)

サブプロットの数をサンプリングごとに確認する
bt[,":="(num_subplot=length(unique(subplot))),by=.(MAIN_PLOT_ID,FinY,sampling)]

グループごとに計算して、新しい変数を作る（congregate subplot abun into the main plot abun）
newbt<-bt[,lapply(.SD,sum),by=.(MAIN_PLOT_ID,FinY,sampling),.SDcols=8:463]​​

列で計算し直す
total_abun<-bt[,sum(.SD),.SDcols=8:463,by=.(MAIN_PLOT_ID,FinY,sampling)]​

nlme::lmeで複数のランダム効果を組み込む（piecewiseSEMなど、lmerが使いにくいとき）
mod1<-lme(Spnbsp_corr~sqrttotal_BA_corr+MATave+CMIave+logCNrt+sqrtSA+Y,
            random=list(PLOT_ID=~1,Prov=~Y),
            control=lmeControl(opt="optim"),can)
summary(mod1)$tT

lmerの警告："Model failed to converge" warning in lmer()
What is "restricted maximum likelihood" and when should it be used?

プロットごとに変数のバリエーションの確認
data[,length(unique(Reg)),by=Plot][,table(V1)] もしくはsummary

グループごとにばらつきを確認する
tapply(can$SA1, can$Prov, summary)
can[, as.list(summary(SA1)), by = Prov]

スペースで区切られた文字列の前部分だけを取り出す（パッケージ stringr​）
CAN[,":="(Genus=str_split(Species,pattern = " ",simplify = TRUE)[,1])]

Notes on multiple random effects in lmer vs lme
#Refered to this (also this), but Multiple random intercepts cannot be included in lme, only works with lmer
m1<-lmer(AGBGI~logtotal_BA+lognbsp+logSA+
           (1|Eco)+(1|PLOT_ID),
         weights=wi,
         can)
summary(m1)$co
m2<-lme(AGBGI~logtotal_BA+lognbsp+logSA,
        random=list(Eco=~1, PLOT_ID=~1),
        weights=~1/wi,
        can)#this only produces the same coefs when modelled with only PLOT_ID
summary(m2)$tT
#but Note that "random=list(Eco=~1)" in lme produces indentical coefs with "(1|Eco)" in lmer

#Nested random intercept (plot id nested to eco) works with lme, same coefs with lmer
m1<-lmer(AGBGI~logtotal_BA+lognbsp+logSA+
           (1|Eco/PLOT_ID),
         weights=wi,
         can)
summary(m1)$co
m2<-lme(AGBGI~logtotal_BA+lognbsp+logSA,
        random=~1|Eco/PLOT_ID,
        weights=~1/wi,
        can)
summary(m2)$tT

geom_bar, 正負の棒グラフごとにアスタリスクを付ける
 geom_text(aes(y = logBAI_AI_mean + 0.005*sign(logBAI_AI_mean), label = "***"), 
            position = position_dodge(width = 0.9), 
            size = 8 , angle = 0)+
※ 上記ではカテゴリごとに分けられないので、その場合は：
  annotate("text", x = 1, y =0.018, label = expression(paste("***")), parse=TRUE, size=6)+
  annotate("text", x = 2, y = -0.008, label = expression(paste("*")), parse=TRUE, size=6)+​

ggplotの図の軸ラベル名で乗算記号を表す​ 　(ギリシャ文字上付き等)
GrowT_mb[,":="(Source=gsub(":","%*%",Source))] %>% 
  .[,":="(Source=gsub("FD","FD[or]*FRed",Source))]#アスタリスク*で下付きエスケープ
ggplot(GrowT_cont,aes(x=mean.Grow,xmin=low.Grow,xmax=up.Grow,y=Source,col=Model))+
scale_y_discrete(labels = scales::parse_format())+

ggplot2で、log(y+1)した応答変数をy軸にscale backして図示する
scale_y_continuous(trans=log1p_trans())​

Year-Month-DayからMonthを抜き出す
temp5060[,":="(Mon=format(as.Date(YMD,format="%Y-%m-%d"),"%m"))]#chr​

図が表示されなくなったら
dev.list()
RStudioGD()

---