How to store summary function into a vector then for loop in R?-CodePudding

I have the following dataset. For each month and site, I'm trying to use a particular package (EcoSim) to estimate overlap (RA4Model).

I want to achieve two things:

install.packages("EcoSimR")
library(EcoSimR)
set.seed(111)
month <- rep(c("J","J","J","F"), each = 4)
site <- rep(c("1","2","3","1"), each = 4)
species <- rep(c("A","B","C","D"), rep = 4)
q1 <- rnorm(16,5,1)
q2 <- rnorm(16,5,1)
q3 <- rnorm(16,5,1)
q4 <- rnorm(16,5,1)
q5 <- rnorm(16,5,1)

df <- data.frame(month, site, species,q1,q2,q3,q4,q5)
df.site <- df[df$month == "J" & df$site == "1",]
df.site <- df.site[,-c(1,2,3)]

RA4model <- niche_null_model(speciesData=df.site, 
                             algo="ra4", metric="pianka", 
                             suppressProg=TRUE,nReps=5000)

summary(RA4model)

First, I want to store certain values from the summary into the corresponding rows in R Here is the output of summary(RA4Model) for Month "J" and site "1"

# Time Stamp:  Sun Jan  8 11:43:05 2023 
# Reproducible:  FALSE 
# Number of Replications:  5000 
# Elapsed Time:  1 secs 
# Metric:  pianka 
# Algorithm:  ra4 
# Observed Index:  0.96516 
# Mean Of Simulated Index:  0.95101 
# Variance Of Simulated Index:  9.4906e-05 
# Lower 95% (1-tail):  0.93823 
# Upper 95% (1-tail):  0.96979 
# Lower 95% (2-tail):  0.93729 
# Upper 95% (2-tail):  0.97272 
# Lower-tail P =  0.8982 
# Upper-tail P =  0.1018 
# Observed metric > 4491 simulated metrics 
# Observed metric < 509 simulated metrics 
# Observed metric = 0 simulated metrics 
# Standardized Effect Size (SES):  1.453

Storing this output into vector. I don't know how to store the lower-1tailP and SES from the summary output into the dataset

df.out <- df[,c(1,2,3)]
df.out$Obs <- RA4model$Obs
df.out$Sim <- mean(RA4model$Sim)
df.out$lower-1tailP <-  #This should be 0.93
df.out$SES <-  #This should be 1.453

Next, I want to loop this so that it does it for each unique(month,site). So the final dataframe looks something like this:

month site Obs Sim lower-1tailP SES
  J     1   ..  ..   ..          ..
  J     2   ..   ..  ..          ..
  J     3   ..   ..  ..          ..
  F     1   ..   ..  ..          ..

CodePudding user response：

You can create a helper function get_eco_sim_result(), which returns a list of the parameters of interest:

get_eco_sim_result <- function(spd, algo= "ra4", metric = "pianka", nReps=500) {
  model = niche_null_model(speciesData = spd,
    algo = algo,metric =metric, nReps = nReps, suppressProg = TRUE
  )
  return(list(
    Obs = model$Obs,
    Sim = mean(model$Sim),
    lower_1tailp = quantile(model$Sim,0.05),
    SES = (model$Obs - mean(model$Sim))/sd(model$Sim)
  ))
}

Then use lapply() to apply that helper function to each subset of the data; here I obtain the subsets of the data using split(). By wrapping the retuned list in data.frame(), you can subsequently use do.call() with rbind()

do.call(
  rbind, lapply(split(df, list(month,site), drop=T), \(d) {
    data.frame(get_eco_sim_result(d[,-c(1,2,3)], nReps=5000))
  })
)

Output:

          Obs       Sim lower_1tailp       SES
F.1 0.9641760 0.9546306    0.9429722 1.0966176
J.1 0.9651613 0.9508969    0.9381335 1.4635265
J.2 0.9931026 0.9842322    0.9800247 2.9524328
J.3 0.9726413 0.9674799    0.9586858 0.7669562

Your original question, however, is perhaps less about the helper function and the loop, but rather about how the parameters in summary(RA4model) are estimated? Use getAnywhere(summary.nichenullmod) to see how these are estimated.