# Nature同种负密度制约论文中的统计错误

Nature同种负密度制约论文中的统计错误

为了弄清这一问题，我重新拟合了其模型，并根据其结果，制作不同水分条件下的

Lebrija-Trejos, E., A. Hernández, and S. J. Wright. 2023. Effects of moisture and density-dependent interactions on tropical tree diversity. Nature 615:100-104.

## LOAD PACKAGES AND FILES ####
library(nlme)
library(lme4)
library(arm)
library(DHARMa)
library(foreach)
library(doParallel)
library(reshape2)
library(vegan)
library(cocor)
library(lubridate)
library(doBy)

bestmodel <- glmer(ALIVE ~ CSD + HSD.fort  + CLD + HLD + SMw  + CSD:SMw + HSD.fort:SMw + SpMO + SpMO:SMw
+ (1 + CSD + CLD + HLD + SMw|SPECIES)
+ (0 + HSD.fort + CSD:SMw + HSD.fort:SMw||SPECIES)
+ (1|PLOT) + (1|TRAP)
+ (1 + CSD + HSD.fort + CLD + HLD +SpMO|YEAR),
control = glmerControl(optCtrl = list(maxfun = 1e+06)),
data = Data, family = binomial, verbose = 0, nAGQ = 1)

############### fig.3 b

SMw<-sort(unique(Data$SMw)) CSD<-seq(0,25,0.1) srv.p.sim<-list() for (i in 1:20){ srv.p.sim[[i]]<-invlogit(fixef(bestmodel)[1] + fixef(bestmodel)[2]*CSD+ fixef(bestmodel)[3]* SMw[i]+ fixef(bestmodel)[4]* SMw[i]*CSD )} mean_p<-invlogit(fixef(bestmodel)[1] + fixef(bestmodel)[2]*CSD+ fixef(bestmodel)[3]* median(SMw)+ fixef(bestmodel)[4]* median(SMw)*CSD) d1<-data.frame(y=srv.p.sim) colnames(d1)<-c(1:20) d1$CSD<-CSD
d2<-gather(d1, key="SMw", value="Pred", 1:20)
d2$moistlevel<-as.numeric(d2$SMw)

pal <- colorRampPalette(c("red","green","blue"))
my_colors <- pal(20)

########## original scale

f3b_v1<-ggplot()+
geom_line(data=d2,aes(x=CSD, y=Pred,group=SMw,color=moistlevel))+
geom_line(aes(x=CSD,y=mean_p), size=1, color="black",
linetype = "dashed")+
theme(legend.position="none")+
ylab("First year survival probability") +
xlab(xlab)+
theme_bw()+
theme(axis.text=element_text(size=12,face="bold"),
axis.title=element_text(size=12,face="bold"))+
theme(strip.text.x = element_text(size =12),
strip.text.y = element_text(size =12),
legend.position='none')

########## logit scale

srv.p.sim2<-list()
for (i in 1:20){
srv.p.sim2[[i]]<-fixef(bestmodel)[1] +
fixef(bestmodel)[2]*CSD+
fixef(bestmodel)[3]* SMw[i]+
fixef(bestmodel)[4]* SMw[i]*CSD
}

d3<-data.frame(y=srv.p.sim2)
colnames(d3)<-c(1:20)
d3$CSD<-CSD d4<-gather(d3, key="SMw", value="Pred_logit_scale", 1:20) d4$moistlevel<-as.numeric(d4\$SMw)

logit_mean_p<-fixef(bestmodel)[1] +
fixef(bestmodel)[2]*CSD+
fixef(bestmodel)[3]* median(SMw)+
fixef(bestmodel)[4]* median(SMw)*CSD

f3b_v2<-ggplot()+
geom_line(data=d4, aes(x=CSD, y=Pred_logit_scale,group=SMw,color=moistlevel))+
geom_line(aes(x=CSD,y=logit_mean_p), size=1, color="black",
linetype = "dashed")+
theme(legend.position="none")+
ylab("logit (First year survival probability)") +
xlab(xlab)+
theme_bw()+
theme(axis.text=element_text(size=12,face="bold"),
axis.title=element_text(size=12,face="bold"))+
theme(strip.text.x = element_text(size =12),
strip.text.y = element_text(size =12),
legend.position='none')

########## logit scale, but relabel the y axis

f3b_v3<-ggplot()+
geom_line(data=d4, aes(x=CSD, y=Pred_logit_scale, group=SMw,color=moistlevel))+
geom_line(aes(x=CSD,y=logit_mean_p), size=1, color="black",
linetype = "dashed")+
scale_y_continuous(breaks = c( logit(0.0001), logit(0.001), logit(0.01),logit(0.1), logit(0.2), logit(0.4), logit(0.6), logit(0.8)),
labels=c(0.0001, 0.001, 0.01, 0.1, 0.2, 0.4, 0.6, 0.8))+
theme(legend.position="none")+
ylab("First year survival probability") +
xlab(xlab)+
theme_bw()+
theme(axis.text=element_text(size=12,face="bold"),
axis.title=element_text(size=12,face="bold"))+
theme(strip.text.x = element_text(size =12),
strip.text.y = element_text(size =12),
legend.position='none')

library(cowplot)
Figure_3b_bestmodel<-plot_grid(f3b_v1, f3b_v2,f3b_v3, nrow=1, ncol=3)
ggsave("Figure_3b_bestmodel.tiff", width=11.61, height=3.42, dpi=600)

https://blog.sciencenet.cn/blog-3442043-1407938.html

## 全部精选博文导读

GMT+8, 2024-2-29 23:36