Algumas questões para seu projeto de pesquisa:
library(optimx)
library(FD)
library(snow)
library(parallel)
library(BioGeoBEARS)
# Podemos fazer um rápido ensaio...
#######################################################
# Statistics -- BAYAREALIKE vs. BAYAREALIKE+J
#######################################################
# We have to extract the log-likelihood differently, depending on the
# version of optim/optimx
LnL_2 = get_LnL_from_BioGeoBEARS_results_object(resBAYAREALIKE)
LnL_1 = get_LnL_from_BioGeoBEARS_results_object(resBAYAREALIKEj)
numparams1 = 3
numparams2 = 2
stats = AICstats_2models(LnL_1, LnL_2, numparams1, numparams2)
stats
# BAYAREALIKE, null model for Likelihood Ratio Test (LRT)
res2 = extract_params_from_BioGeoBEARS_results_object(results_object=resBAYAREALIKE, returnwhat="table", addl_params=c("j"), paramsstr_digits=4)
# BAYAREALIKE+J, alternative model for Likelihood Ratio Test (LRT)
res1 = extract_params_from_BioGeoBEARS_results_object(results_object=resBAYAREALIKEj, returnwhat="table", addl_params=c("j"), paramsstr_digits=4)
rbind(res2, res1)
conditional_format_table(stats)
tmp_tests = conditional_format_table(stats)
restable = rbind(restable, res2, res1)
teststable = rbind(teststable, tmp_tests)
#####################################################################################
#####################################################################################
# carregando os dados:
res = "Psychotria_DIVALIKE_M0_unconstrained_v1.Rdata"
load(res)
res_onlyGeo_divaLike = res
res = "Psychotria_DIVALIKE+J_M0_unconstrained_v1.Rdata"
load(res)
res_onlyGeo_divaLikej = res
res = "Psychotria_BAYAREALIKE_M0_unconstrained_v1.Rdata"
load(res)
res_onlyGeo_bayareaLike = res
res = "Psychotria_BAYAREALIKE+J_M0_unconstrained_v1.Rdata"
load(res)
res_onlyGeo_bayareaLikej = res
res = "Psychotria_DEC_M0_unconstrained_v1.Rdata"
load(res)
res_onlyGeo_dec = res
res = "Psychotria_DEC+J_M0_unconstrained_v1.Rdata"
load(res)
res_onlyGeo_decj = res
#######################################
# para calcular o log da verossimilhança, é necessário o número de parâmetros de cada modelo:
numParam_Bay_M3 <- 2
numParam_Bayj_M3 <- 3
numParam_Diva_M3 <- 2
numParam_Divaj_M3 <- 3
numParam_DEC_M3 <- 2
numParam_DECj_M3 <- 3
#######################################
# lista de resultados:
results <- list()
results[[1]] <- res_onlyGeo_divaLike; results[[2]] <- res_onlyGeo_divaLikej
results[[3]] <- res_onlyGeo_bayareaLike; results[[4]] <- res_onlyGeo_bayareaLikej
results[[5]] <- res_onlyGeo_dec; results[[6]] <- res_onlyGeo_decj
# Set up empty tables to hold the statistical results
restable = NULL
teststable = NULL
resPrel <- list()
for(i in 1:length(results)){
resPrel[[i]] = extract_params_from_BioGeoBEARS_results_object(results_object = results[[i]],
returnwhat = "table", addl_params = c("j"), paramsstr_digits = 4)
}
resPrel
# modelo 1 e modelo 2 são aninhados:
stats = AICstats_2models(resPrel[[2]]$LnL, resPrel[[1]]$LnL, 5, 4)
stats
tmp_tests = conditional_format_table(stats)
tmp_tests
restable = rbind(restable, resPrel[[2]]$LnL, resPrel[[1]]$LnL)
teststable = rbind(teststable, tmp_tests)
teststable
# modelo 3 e modelo 4:
stats = AICstats_2models(resPrel[[4]]$LnL, resPrel[[3]]$LnL, 5, 4)
stats
tmp_tests = conditional_format_table(stats)
tmp_tests
restable = rbind(restable, resPrel[[4]]$LnL, resPrel[[3]]$LnL)
teststable = rbind(teststable, tmp_tests)
teststable
# modelo 5 e modelo 6:
stats = AICstats_2models(resPrel[[6]]$LnL, resPrel[[5]]$LnL, 5, 4)
stats
tmp_tests = conditional_format_table(stats)
tmp_tests
restable = rbind(restable, resPrel[[6]]$LnL, resPrel[[5]]$LnL)
teststable = rbind(teststable, tmp_tests)
teststable
######################################################################
#########################################################################
# ASSEMBLE RESULTS TABLES: DEC, DEC+J, DIVALIKE, DIVALIKE+J, BAYAREALIKE, BAYAREALIKE+J
#########################################################################
teststable$alt = c("DIVALIKE+J", "BAYAREALIKE+J", "DEC+J")
teststable
teststable$null = c("DIVALIKE", "BAYAREALIKE", "DEC")
teststable
dim(restable)
row.names(restable) = c("DIVALIKE+j", "DIVALIKE", "BAYAREALIKE+j", "BAYAREALIKE",
"DEC+j", "DEC")
restable = put_jcol_after_ecol(restable)
restable = round(restable, 2)
colnames(restable) <- 'lnL'
# Look at the results!!
restable
teststable
Algumas questões para seu projeto de pesquisa:
#######################################################
# Pesos dos modelos
#######################################################
restable
dim(restable)
# restable2 = as.matrix(restable[-c(1:12), ])
restable2 = as.matrix(restable)
restable2
dim(restable2)
teststable
# With AICs:
?calc_AIC_column
modelos <- as.numeric(unlist(c(teststable$LnLalt[1:3], teststable$LnLnull[1:3])))
modelos
# ordem do restable2:
modelos2 <- as.numeric(unlist(c(teststable$LnLalt[1], teststable$LnLnull[1],
teststable$LnLalt[2], teststable$LnLnull[2],
teststable$LnLalt[3], teststable$LnLnull[3])))
modelos2
# mesma ordem do restable2
gl2 <- as.numeric(unlist(c(teststable$DFalt[1], teststable$DFnull[1],
teststable$DFalt[2], teststable$DFnull[2],
teststable$DFalt[3], teststable$DFnull[3])))
gl2 # na mesma ordem do restable2
AICtable2 = calc_AIC_column(LnL_vals = modelos2,
nparam_vals = gl2)
restable_soma <- modelos2
restable_soma
restable_soma = as.matrix(restable_soma)
colnames(restable_soma) <- 'lnL'
row.names(restable_soma) = c("DIVALIKE+j", "DIVALIKE", "BAYAREALIKE+j", "BAYAREALIKE",
"DEC+j", "DEC")
restable_soma
restable_soma = cbind(restable_soma, AICtable2)
restable_soma
restable_AIC_rellike = AkaikeWeights_on_summary_table(restable=restable_soma, colname_to_use="AIC")
restable_AIC_rellike = put_jcol_after_ecol(restable_AIC_rellike)
conditional_format_table(restable_AIC_rellike)
# se eu quiser corrigir porque tenho poucas espécies:
# With AICcs -- factors in sample size
samplesize = 38
AICtable = round(calc_AICc_column(modelos2,
nparam_vals = gl2, samplesize = samplesize), 4)
AICtable
restable2 = cbind(round(restable_soma, 2), AICtable)
restable2
restable_AICc_rellike = AkaikeWeights_on_summary_table(restable=restable2, colname_to_use="AICc")
restable_AICc_rellike = put_jcol_after_ecol(restable_AICc_rellike)
conditional_format_table(restable_AICc_rellike)
restable_AICc_rellike <- round(restable_AICc_rellike, 4)
restable_AICc_rellike