Hi all,
Recently, I’m doing some analysis about the influence of biotic neighborhood on the survival of seedling. I have seedling survival data from 450 1×1m plot across 11 years.
The main question I want to resolve is:
- how different neighbor crowding( e.g. number of conspecific seedling, number of heterospecific seedling, the sum of basal area of conspecific adult, the sum of basal area of heterospecific adult) influence seedling survival?
- how functional trait can mediate species response different neighbor crowding effect?
When I ran the data based on the hierarchical Bayesian model, none of the above factors can have significant influence( with 95% intervals don’t overlap with 0). As people in my study area have found those factor do have influence, so I try the generalized linear mixed model (GLMM) based on lme4 package and stan_glmer() based on rstanrm package. Both methods give me a significant result.
So which one should I trust? I don’t know what happens. I’m wondering if anyone can give me a hand. Thank you so much!
Following is the hierarchical Bayesian model and two other methods I have tried.
-
tree.all.yr.c.sp = glmer(survival ~ (1|year) + (1|code) + Shom + Shet + Acon + Ahet, family = binomial, data= tree.dat)
-
rstanarm::stan_glmer(Status ~ (1|plot) + (1|census) + Shom + Shet + Acon + Ahet,
data = seedling.comp, family = binomial, QR = TRUE,
# this next line is only to keep the example small in size!
chains = 3, cores = 3, seed = 12345, iter = 500)
data {
int<lower=0> nplot; // number of seedling plot
int<lower=0> nspe; // number of species
int<lower=0> ncensus;//number of census time
int<lower=0> n; //number of survival record
int<lower=0> ncov; //number of covariates
int<lower=0> plot[n]; //plot name
int<lower=0> spe[n];//species name
int<lower=0> census[n];//census time, year
vector[nspe] meanT; //functional trait
matrix[n,ncov] COV; // 1. intercept; 2.conspecific seedling; 3.hetspecific seedling; 4.conspecific adult; 5. hetspecific adult
int<lower=0,upper=1> status[n]; //status of each record
}
parameters {
real bet[ncov];
real b[nspe,ncov];
real int1[ncov];
real int2[nplot];
real int3[ncensus];
real <lower=0.00001> sigmab[5];
real <lower=0.00001> sigmaplot;
real <lower=0.00001> sigmacensus;
}
model {
sigmacensus~gamma(100,100);
sigmaplot~gamma(100,100);
for(i in 1:5)
{
bet[i]~normal(0,100);
int1[i]~normal(0,100);
sigmab[i]~gamma(100,100);
}
for(i in 1:nplot)
{
int2[i]~normal(0,sigmaplot);
}
for(i in 1:ncensus)
{
int3[i]~normal(0,sigmacensus);
}
for(j in 1:nspe)
{
b[j,1]~normal(int1[1]+bet[1]*meanT[j],sigmab[1]); //bet[1] represent how functional trait related to seedling survival
b[j,2]~normal(int1[2]+bet[2]*meanT[j],sigmab[2]); //bet[2] represent how trait mediated the effect of conspecific seedling
b[j,3]~normal(int1[3]+bet[3]*meanT[j],sigmab[3]); //bet[2] represent how trait mediated the effect of hetspecific seedling
b[j,4]~normal(int1[4]+bet[4]*meanT[j],sigmab[4]); //bet[2] represent how trait mediated the effect of conspecific adult tree
b[j,5]~normal(int1[5]+bet[5]*meanT[j],sigmab[5]); //bet[2] represent how trait mediated the effect of hetspecific adult tree
}
for(i in 1:n)
status[i]~bernoulli_logit(b[spe[i],1]*COV[i,1]+b[spe[i],2]*COV[i,2]+b[spe[i],3]*COV[i,3]+b[spe[i],4]*COV[i,4]+b[spe[i],5]*COV[i,5]+int2[plot[i]]+int3[census[i]]);
}
****