Messika et al. 2017: R scripts for stage 1 of the statistical analyses
Here are presented the R scripts used in Messika et al.’s 2017 project on determining factors influencing the reproductive success of arthropod vectors on host-parasite communities in the western Negev Desert in Israel.
Briefly, data analyses were conducted in two stages. First, we searched for the most important factors that best explained the variability in female quality and reproductive success (Table 2), and then, we quantified the causal pathways of the intrinsic and extrinsic factors, using a path analysis approach (Wootton 1994).
In stage 1, we explored 13 model-sets, each with a different dependent variable. In 10 of the model-sets, the dependent variable is related to the current reproductive success of females, while in three of the model-sets, it is related to female quality. Each model-set included all possible combinations of additive nested GLM models. For example, if a model-set included three variables (V1, V2 and V3), we fitted a total of eight GLMs including: null model, V1, V2, V3, V1+V2, V1+V3, V2+V3, V1+V2+V3. We scored each model using AICc and calculated their AICc weight (wAICc, Burnham and Anderson 2002). Then, we assessed the direction of explanatory variables based on the averaged regression coefficients, using wAICc as weights. As the model-set included all possible combinations of GLMs (excluding interactions between variables), each candidate predictor variable appears in the same number of models (Burnham and Anderson 2002, p. 169). Therefore, the sum of the wAICc of all models that included a certain focal variable can be used as a measure of relative variable importance (Burnham and Anderson 2002; Johnson and Omland 2004). In all analyses at this stage, we used the function ‘glm’ in the ‘AICcmodavg’ package in R (R-Core-Team 2013).
Beloware presented the R scripts used to calculate the AICc of different models and the relative importance of each variable.
Data file used in present R scripts: Messika et al.’s raw data datasheet.xlsx (archived in Figshare.com:
To run the scripts, data must be transformed to .csv (comma delimited) file.
R scripts
1. Determinants of the variability of female fleas´ offspring number
####################################################
#####################################################
# clear prev objects
rm(list=ls())
#####################################################
# set working directory
#setwd("~/Work/tree of knowledge/1-manuscripts/4- in preperation/Jost partioning/analysis/GLMs")
#####################################################
# library packages
prep.packages = function(package.list)
{
loaded = package.list %in% .packages()
if ( all(loaded) ) return(invisible())
package.list = package.list[!loaded]
installed = package.list %in% .packages(TRUE)
if ( !all(installed) ) install.packages(package.list[!installed])
for ( p in package.list )
{
print(paste("Loading package:",p))
library(p,character.only=TRUE)
}
}
prep.packages(c("AICcmodavg","gdata"))
#####################################################
# name of the general input data file
Input_1_name= file.choose("Data") # name of the input data table as text file
# set working directory
setwd(dirname(Input_1_name))
# names of output files
# The name of the output AIC table
Output_name_1 = paste("AIC_Tab.txt")
# the name of the model list output table
Output_name_2 = paste("model_list.txt")
# the name of the importance and model averaging table
Output_name_3 = paste("Impo_mod_ave.txt")
#####################################################
#####################################################
# read the data
Main_Data = read.csv(Input_1_name)
######################################################
# here is where you can subset the main data
# if only a portion of the data points or variables are needed
#here is an example:
# work_data=subset(work_data,work_data$Cont_in==mainland_in) # subset according to mainliand_in
work_data=Main_Data
names_var=names(work_data) # names of variables
#### immportant #####
# Change the numbers within the [] to the colimns that correspond to exlanatory and reponse variables
names_expln=names_var[c(21,23,38:43,50:52)] # names of explanatory variables
names_res=names_var[5] # names of response variables
#### change here below according to the number of exlanatory variables in your dataset -
# change the 5 to the ocrrect value
Number_Of_Explan= 11
for ( i in 38:39 ) work_data[[i]] = as.factor(work_data[[i]])
for ( i in 51:52 ) work_data[[i]] = as.ordered(work_data[[i]])
for ( i in 50 ) work_data[[i]] = as.ordered(work_data[[i]])
#####################################################
#### important - change her to focus on the reponse variable you want to work on
res_var=1 # The response variable you want to work on now
#####################################################
# create the model list
n <- length(names_expln)
id <- unlist(
lapply(1:n,
function(i)combn(1:n,i,simplify=F)
)
,recursive=F)
gen.form_1<-function(x) as.formula(paste(paste(names_res[res_var]),'~ ',paste(names_expln[x],collapse='+')))
formulas_alpha<-lapply(id, gen.form_1)
num_models=length(id)
Cand.mod <- list()
for (j in 1:num_models)
{
Cand.mod[[j]]=glm(formulas_alpha[[j]],family= "poisson", data=work_data)
}
#####################################################
# add the null glm - no variables at all - needed for the importance to balance the design
formula_null=as.formula(paste(paste(names_res[res_var]),'~ 1'))
Cand.mod[[(num_models+1)]]=glm(formula_null,family="poisson", data=work_data)
#####################################################
##create a vector of names to trace back models in set
Modnames <- paste("mod", 1:length(Cand.mod), sep = "_")
#####################################################
## join the null with other models
formulas_alpha_1=c(formulas_alpha,formula_null)
#####################################################
# the AIcc table
results_1=aictab(cand.set = Cand.mod, modnames = Modnames, sort = FALSE)
#####################################################
# the model list table
results_2=cbind(Modnames,formulas_alpha_1)
#####################################################
names(Cand.mod)<-Modnames
attach(Cand.mod)
results_3 = array(0,c(Number_Of_Explan,6))
results_3 = as.data.frame(results_3)
names(results_3)=list("Var","Importance","Coef_est","Coef_SE","Coef_Low_Cf","Coef_high_cf")
names_coef<-names(coef(mod_2047))
names_coef<-names_coef[names_coef!="(Intercept)"]
for(count_1 in 1:Number_Of_Explan)
{
Imp_1=importance(cand.set = Cand.mod, parm = names_coef[count_1], modnames = Modnames,
second.ord = TRUE, nobs = NULL)
results_3[count_1,1]=Imp_1$parm
results_3[count_1,2]=Imp_1$w.plus
mod_ave_1=modavg(parm = names_coef[count_1], cand.set = Cand.mod, modnames = Modnames)
results_3[count_1,3]=mod_ave_1$Mod.avg.beta
results_3[count_1,4]=mod_ave_1$Uncond.SE
results_3[count_1,5]=mod_ave_1$Lower.CL
results_3[count_1,6]=mod_ave_1$Upper.CL
}
#####################################################
#write the output tables
write.table(results_1,Output_name_1, row.names=F)
write.table(results_2,Output_name_2, row.names=F)
write.table(results_3,Output_name_3, row.names=F)
#write.table(results_2,"model_list_1.txt", row.names=F,sep=";")
#####################################################
# end of script
#####################################################
2. Determinants of the variability of the residual values of the regression of offspring number on mother mass
####################################################
#####################################################
# clear prev objects
rm(list=ls())
#####################################################
# set working directory
#setwd("~/Work/tree of knowledge/1-manuscripts/4- in preperation/Jost partioning/analysis/GLMs")
#####################################################
# library packages
prep.packages = function(package.list)
{
loaded = package.list %in% .packages()
if ( all(loaded) ) return(invisible())
package.list = package.list[!loaded]
installed = package.list %in% .packages(TRUE)
if ( !all(installed) ) install.packages(package.list[!installed])
for ( p in package.list )
{
print(paste("Loading package:",p))
library(p,character.only=TRUE)
}
}
prep.packages(c("AICcmodavg","gdata"))
#####################################################
# name of the general input data file
Input_1_name= file.choose("Data") # name of the input data table as text file
# set working directory
setwd(dirname(Input_1_name))
# names of output files
# The name of the output AIC table
Output_name_1 = paste("AIC_Tab.txt")
# the name of the model list output table
Output_name_2 = paste("model_list.txt")
# the name of the importance and model averaging table
Output_name_3 = paste("Impo_mod_ave.txt")
#####################################################
#####################################################
# read the data
Main_Data = read.csv(Input_1_name)
######################################################
# here is where you can subset the main data
# if only a portion of the data points or variables are needed
#here is an example:
# work_data=subset(work_data,work_data$Cont_in==mainland_in) # subset according to mainliand_in
work_data=Main_Data
names_var=names(work_data) # names of variables
#### immportant #####
# Change the numbers within the [] to the colimns that correspond to exlanatory and reponse variables
names_expln=names_var[c(50,43,38,39,41,42,40,23,51,52)] # names of explanatory variables
names_res=names_var[53] # names of response variables
#### change here below according to the number of exlanatory variables in your dataset -
# change the 5 to the ocrrect value
Number_Of_Explan= 10
for ( i in 38:39 ) work_data[[i]] = as.factor(work_data[[i]])
for ( i in 51:52 ) work_data[[i]] = as.ordered(work_data[[i]])
for ( i in 50 ) work_data[[i]] = as.ordered(work_data[[i]])
#####################################################
#### important - change her to focus on the reponse variable you want to work on
res_var=1 # The response variable you want to work on now
#####################################################
# create the model list
n <- length(names_expln)
id <- unlist(
lapply(1:n,
function(i)combn(1:n,i,simplify=F)
)
,recursive=F)
gen.form_1<-function(x) as.formula(paste(paste(names_res[res_var]),'~ ',paste(names_expln[x],collapse='+')))
formulas_alpha<-lapply(id, gen.form_1)
num_models=length(id)
Cand.mod <- list()
for (j in 1:num_models)
{
Cand.mod[[j]]=glm(formulas_alpha[[j]],family= "gaussian", data=work_data)
}
#####################################################
# add the null glm - no variables at all - needed for the importance to balance the design
formula_null=as.formula(paste(paste(names_res[res_var]),'~ 1'))
Cand.mod[[(num_models+1)]]=glm(formula_null,family="gaussian", data=work_data)
#####################################################
##create a vector of names to trace back models in set
Modnames <- paste("mod", 1:length(Cand.mod), sep = "_")
#####################################################
## join the null with other models
formulas_alpha_1=c(formulas_alpha,formula_null)
#####################################################
# the AIcc table
results_1=aictab(cand.set = Cand.mod, modnames = Modnames, sort = FALSE)
#####################################################
# the model list table
results_2=cbind(Modnames,formulas_alpha_1)
#####################################################
names(Cand.mod)<-Modnames
attach(Cand.mod)
results_3 = array(0,c(Number_Of_Explan,6))
results_3 = as.data.frame(results_3)
names(results_3)=list("Var","Importance","Coef_est","Coef_SE","Coef_Low_Cf","Coef_high_cf")
names_coef<-names(coef(mod_1023))
names_coef<-names_coef[names_coef!="(Intercept)"]
for(count_1 in 1:Number_Of_Explan)
{
Imp_1=importance(cand.set = Cand.mod, parm = names_coef[count_1], modnames = Modnames,
second.ord = TRUE, nobs = NULL)
results_3[count_1,1]=Imp_1$parm
results_3[count_1,2]=Imp_1$w.plus
mod_ave_1=modavg(parm = names_coef[count_1], cand.set = Cand.mod, modnames = Modnames)
results_3[count_1,3]=mod_ave_1$Mod.avg.beta
results_3[count_1,4]=mod_ave_1$Uncond.SE
results_3[count_1,5]=mod_ave_1$Lower.CL
results_3[count_1,6]=mod_ave_1$Upper.CL
}
#####################################################
#write the output tables
write.table(results_1,Output_name_1, row.names=F)
write.table(results_2,Output_name_2, row.names=F)
write.table(results_3,Output_name_3, row.names=F)
#write.table(results_2,"model_list_1.txt", row.names=F,sep=";")
#####################################################
# end of script
#####################################################
3. Determinants of the variability of offspring sex ratio
####################################################
#####################################################
# clear prev objects
rm(list=ls())
#####################################################
# set working directory
#setwd("~/Work/tree of knowledge/1-manuscripts/4- in preperation/Jost partioning/analysis/GLMs")
#####################################################
# library packages
prep.packages = function(package.list)
{
loaded = package.list %in% .packages()
if ( all(loaded) ) return(invisible())
package.list = package.list[!loaded]
installed = package.list %in% .packages(TRUE)
if ( !all(installed) ) install.packages(package.list[!installed])
for ( p in package.list )
{
print(paste("Loading package:",p))
library(p,character.only=TRUE)
}
}
prep.packages(c("AICcmodavg","gdata"))
#####################################################
# name of the general input data file
Input_1_name= file.choose("Data") # name of the input data table as text file
# set working directory
setwd(dirname(Input_1_name))
# names of output files
# The name of the output AIC table
Output_name_1 = paste("AIC_Tab.txt")
# the name of the model list output table
Output_name_2 = paste("model_list.txt")
# the name of the importance and model averaging table
Output_name_3 = paste("Impo_mod_ave.txt")
#####################################################
#####################################################
# read the data
Main_Data = read.csv(Input_1_name)
######################################################
# here is where you can subset the main data
# if only a portion of the data points or variables are needed
#here is an example:
# work_data=subset(work_data,work_data$Cont_in==mainland_in) # subset according to mainliand_in
work_data=Main_Data
names_var=names(work_data) # names of variables
#### immportant #####
# Change the numbers within the [] to the colimns that correspond to exlanatory and reponse variables
names_expln=names_var[c(21,23,38:43,50:52)] # names of explanatory variables
names_res=names_var[29] # names of response variables
#### change here below according to the number of exlanatory variables in your dataset -
# change the 5 to the ocrrect value
Number_Of_Explan= 11
for ( i in 29 ) work_data[[i]] = as.factor(work_data[[i]])
for ( i in 38:39 ) work_data[[i]] = as.factor(work_data[[i]])
for ( i in 51:52 ) work_data[[i]] = as.ordered(work_data[[i]])
for ( i in 50 ) work_data[[i]] = as.ordered(work_data[[i]])
#####################################################
#### important - change her to focus on the reponse variable you want to work on
res_var=1 # The response variable you want to work on now
#####################################################
# create the model list
n <- length(names_expln)
id <- unlist(
lapply(1:n,
function(i)combn(1:n,i,simplify=F)
)
,recursive=F)
gen.form_1<-function(x) as.formula(paste(paste(names_res[res_var]),'~ ',paste(names_expln[x],collapse='+')))
formulas_alpha<-lapply(id, gen.form_1)
num_models=length(id)
Cand.mod <- list()
for (j in 1:num_models)
{
Cand.mod[[j]]=glm(formulas_alpha[[j]],family= "binomial", data=work_data)
}
#####################################################
# add the null glm - no variables at all - needed for the importance to balance the design
formula_null=as.formula(paste(paste(names_res[res_var]),'~ 1'))
Cand.mod[[(num_models+1)]]=glm(formula_null,family="binomial", data=work_data)
#####################################################
##create a vector of names to trace back models in set
Modnames <- paste("mod", 1:length(Cand.mod), sep = "_")
#####################################################
## join the null with other models
formulas_alpha_1=c(formulas_alpha,formula_null)
#####################################################
# the AIcc table
results_1=aictab(cand.set = Cand.mod, modnames = Modnames, sort = FALSE)
#####################################################
# the model list table
results_2=cbind(Modnames,formulas_alpha_1)
#####################################################
names(Cand.mod)<-Modnames
attach(Cand.mod)
results_3 = array(0,c(Number_Of_Explan,6))
results_3 = as.data.frame(results_3)
names(results_3)=list("Var","Importance","Coef_est","Coef_SE","Coef_Low_Cf","Coef_high_cf")
names_coef<-names(coef(mod_2047))
names_coef<-names_coef[names_coef!="(Intercept)"]
for(count_1 in 1:Number_Of_Explan)
{
Imp_1=importance(cand.set = Cand.mod, parm = names_coef[count_1], modnames = Modnames,
second.ord = TRUE, nobs = NULL)
results_3[count_1,1]=Imp_1$parm
results_3[count_1,2]=Imp_1$w.plus
mod_ave_1=modavg(parm = names_coef[count_1], cand.set = Cand.mod, modnames = Modnames)
results_3[count_1,3]=mod_ave_1$Mod.avg.beta
results_3[count_1,4]=mod_ave_1$Uncond.SE
results_3[count_1,5]=mod_ave_1$Lower.CL
results_3[count_1,6]=mod_ave_1$Upper.CL
}
#####################################################
#write the output tables
write.table(results_1,Output_name_1, row.names=F)
write.table(results_2,Output_name_2, row.names=F)
write.table(results_3,Output_name_3, row.names=F)
#write.table(results_2,"model_list_1.txt", row.names=F,sep=";")
#####################################################
# end of script
#####################################################
4. Determinants of the variability of male offspring mean development time
####################################################
#####################################################
# clear prev objects
rm(list=ls())
#####################################################
# set working directory
#setwd("~/Work/tree of knowledge/1-manuscripts/4- in preperation/Jost partioning/analysis/GLMs")
#####################################################
# library packages
prep.packages = function(package.list)
{
loaded = package.list %in% .packages()
if ( all(loaded) ) return(invisible())
package.list = package.list[!loaded]
installed = package.list %in% .packages(TRUE)
if ( !all(installed) ) install.packages(package.list[!installed])
for ( p in package.list )
{
print(paste("Loading package:",p))
library(p,character.only=TRUE)
}
}
prep.packages(c("AICcmodavg","gdata"))
#####################################################
# name of the general input data file
Input_1_name= file.choose("Data") # name of the input data table as text file
# set working directory
setwd(dirname(Input_1_name))
# names of output files
# The name of the output AIC table
Output_name_1 = paste("AIC_Tab.txt")
# the name of the model list output table
Output_name_2 = paste("model_list.txt")
# the name of the importance and model averaging table
Output_name_3 = paste("Impo_mod_ave.txt")
#####################################################
#####################################################
# read the data
Main_Data = read.csv(Input_1_name)
######################################################
# here is where you can subset the main data
# if only a portion of the data points or variables are needed
#here is an example:
# work_data=subset(work_data,work_data$Cont_in==mainland_in) # subset according to mainliand_in
work_data=Main_Data
names_var=names(work_data) # names of variables
#### immportant #####
# Change the numbers within the [] to the colimns that correspond to exlanatory and reponse variables
names_expln=names_var[c(5,21,38,43,50,51)] # names of explanatory variables
names_res=names_var[36] # names of response variables
#### change here below according to the number of exlanatory variables in your dataset -
# change the 5 to the ocrrect value
Number_Of_Explan= 6
for ( i in 29 ) work_data[[i]] = as.factor(work_data[[i]])
for ( i in 38:39 ) work_data[[i]] = as.factor(work_data[[i]])
for ( i in 51:52 ) work_data[[i]] = as.ordered(work_data[[i]])
for ( i in 50 ) work_data[[i]] = as.ordered(work_data[[i]])
#####################################################
#### important - change her to focus on the reponse variable you want to work on
res_var=1 # The response variable you want to work on now
#####################################################
# create the model list
n <- length(names_expln)
id <- unlist(
lapply(1:n,
function(i)combn(1:n,i,simplify=F)
)
,recursive=F)
gen.form_1<-function(x) as.formula(paste(paste(names_res[res_var]),'~ ',paste(names_expln[x],collapse='+')))
formulas_alpha<-lapply(id, gen.form_1)
num_models=length(id)
Cand.mod <- list()
for (j in 1:num_models)
{
Cand.mod[[j]]=glm(formulas_alpha[[j]],family= "binomial", data=work_data)
}
#####################################################
# add the null glm - no variables at all - needed for the importance to balance the design
formula_null=as.formula(paste(paste(names_res[res_var]),'~ 1'))
Cand.mod[[(num_models+1)]]=glm(formula_null,family="binomial", data=work_data)
#####################################################
##create a vector of names to trace back models in set
Modnames <- paste("mod", 1:length(Cand.mod), sep = "_")
#####################################################
## join the null with other models
formulas_alpha_1=c(formulas_alpha,formula_null)
#####################################################
# the AIcc table
results_1=aictab(cand.set = Cand.mod, modnames = Modnames, sort = FALSE)
#####################################################
# the model list table
results_2=cbind(Modnames,formulas_alpha_1)
#####################################################
names(Cand.mod)<-Modnames
attach(Cand.mod)
results_3 = array(0,c(Number_Of_Explan,6))
results_3 = as.data.frame(results_3)
names(results_3)=list("Var","Importance","Coef_est","Coef_SE","Coef_Low_Cf","Coef_high_cf")
names_coef<-names(coef(mod_63))
names_coef<-names_coef[names_coef!="(Intercept)"]
for(count_1 in 1:Number_Of_Explan)
{
Imp_1=importance(cand.set = Cand.mod, parm = names_coef[count_1], modnames = Modnames,
second.ord = TRUE, nobs = NULL)
results_3[count_1,1]=Imp_1$parm
results_3[count_1,2]=Imp_1$w.plus
mod_ave_1=modavg(parm = names_coef[count_1], cand.set = Cand.mod, modnames = Modnames)
results_3[count_1,3]=mod_ave_1$Mod.avg.beta
results_3[count_1,4]=mod_ave_1$Uncond.SE
results_3[count_1,5]=mod_ave_1$Lower.CL
results_3[count_1,6]=mod_ave_1$Upper.CL
}
#####################################################
#write the output tables
write.table(results_1,Output_name_1, row.names=F)
write.table(results_2,Output_name_2, row.names=F)
write.table(results_3,Output_name_3, row.names=F)
#write.table(results_2,"model_list_1.txt", row.names=F,sep=";")
#####################################################
# end of script
#####################################################
5. Determinants of the variability of female offspring mean development time
####################################################
#####################################################
# clear prev objects
rm(list=ls())
#####################################################
# set working directory
#setwd("~/Work/tree of knowledge/1-manuscripts/4- in preperation/Jost partioning/analysis/GLMs")
#####################################################
# library packages
prep.packages = function(package.list)
{
loaded = package.list %in% .packages()
if ( all(loaded) ) return(invisible())
package.list = package.list[!loaded]
installed = package.list %in% .packages(TRUE)
if ( !all(installed) ) install.packages(package.list[!installed])
for ( p in package.list )
{
print(paste("Loading package:",p))
library(p,character.only=TRUE)
}
}
prep.packages(c("AICcmodavg","gdata"))
#####################################################
# name of the general input data file
Input_1_name= file.choose("Data") # name of the input data table as text file
# set working directory
setwd(dirname(Input_1_name))
# names of output files
# The name of the output AIC table
Output_name_1 = paste("AIC_Tab.txt")
# the name of the model list output table
Output_name_2 = paste("model_list.txt")
# the name of the importance and model averaging table
Output_name_3 = paste("Impo_mod_ave.txt")
#####################################################
#####################################################
# read the data
Main_Data = read.csv(Input_1_name)
######################################################
# here is where you can subset the main data
# if only a portion of the data points or variables are needed
#here is an example:
# work_data=subset(work_data,work_data$Cont_in==mainland_in) # subset according to mainliand_in
work_data=Main_Data
names_var=names(work_data) # names of variables
#### immportant #####
# Change the numbers within the [] to the colimns that correspond to exlanatory and reponse variables
names_expln=names_var[c(5,21,38,43,50,51)] # names of explanatory variables
names_res=names_var[37] # names of response variables
#### change here below according to the number of exlanatory variables in your dataset -
# change the 5 to the ocrrect value
Number_Of_Explan= 6
for ( i in 29 ) work_data[[i]] = as.factor(work_data[[i]])
for ( i in 38:39 ) work_data[[i]] = as.factor(work_data[[i]])
for ( i in 51:52 ) work_data[[i]] = as.ordered(work_data[[i]])
for ( i in 50 ) work_data[[i]] = as.ordered(work_data[[i]])
#####################################################
#### important - change her to focus on the reponse variable you want to work on
res_var=1 # The response variable you want to work on now
#####################################################
# create the model list
n <- length(names_expln)
id <- unlist(
lapply(1:n,
function(i)combn(1:n,i,simplify=F)
)
,recursive=F)
gen.form_1<-function(x) as.formula(paste(paste(names_res[res_var]),'~ ',paste(names_expln[x],collapse='+')))
formulas_alpha<-lapply(id, gen.form_1)
num_models=length(id)
Cand.mod <- list()
for (j in 1:num_models)
{
Cand.mod[[j]]=glm(formulas_alpha[[j]],family= "binomial", data=work_data)
}
#####################################################
# add the null glm - no variables at all - needed for the importance to balance the design
formula_null=as.formula(paste(paste(names_res[res_var]),'~ 1'))
Cand.mod[[(num_models+1)]]=glm(formula_null,family="binomial", data=work_data)
#####################################################
##create a vector of names to trace back models in set
Modnames <- paste("mod", 1:length(Cand.mod), sep = "_")
#####################################################
## join the null with other models
formulas_alpha_1=c(formulas_alpha,formula_null)
#####################################################
# the AIcc table
results_1=aictab(cand.set = Cand.mod, modnames = Modnames, sort = FALSE)
#####################################################
# the model list table
results_2=cbind(Modnames,formulas_alpha_1)
#####################################################
names(Cand.mod)<-Modnames
attach(Cand.mod)
results_3 = array(0,c(Number_Of_Explan,6))
results_3 = as.data.frame(results_3)
names(results_3)=list("Var","Importance","Coef_est","Coef_SE","Coef_Low_Cf","Coef_high_cf")
names_coef<-names(coef(mod_63))
names_coef<-names_coef[names_coef!="(Intercept)"]
for(count_1 in 1:Number_Of_Explan)
{
Imp_1=importance(cand.set = Cand.mod, parm = names_coef[count_1], modnames = Modnames,
second.ord = TRUE, nobs = NULL)
results_3[count_1,1]=Imp_1$parm
results_3[count_1,2]=Imp_1$w.plus
mod_ave_1=modavg(parm = names_coef[count_1], cand.set = Cand.mod, modnames = Modnames)
results_3[count_1,3]=mod_ave_1$Mod.avg.beta
results_3[count_1,4]=mod_ave_1$Uncond.SE
results_3[count_1,5]=mod_ave_1$Lower.CL
results_3[count_1,6]=mod_ave_1$Upper.CL
}
#####################################################
#write the output tables
write.table(results_1,Output_name_1, row.names=F)
write.table(results_2,Output_name_2, row.names=F)
write.table(results_3,Output_name_3, row.names=F)
#write.table(results_2,"model_list_1.txt", row.names=F,sep=";")
#####################################################
# end of script
#####################################################
6. Determinants of the variability of male offspring average mass
####################################################
#####################################################
# clear prev objects
rm(list=ls())
#####################################################
# set working directory
#setwd("~/Work/tree of knowledge/1-manuscripts/4- in preperation/Jost partioning/analysis/GLMs")
#####################################################
# library packages
prep.packages = function(package.list)
{
loaded = package.list %in% .packages()
if ( all(loaded) ) return(invisible())
package.list = package.list[!loaded]
installed = package.list %in% .packages(TRUE)
if ( !all(installed) ) install.packages(package.list[!installed])
for ( p in package.list )
{
print(paste("Loading package:",p))
library(p,character.only=TRUE)
}
}
prep.packages(c("AICcmodavg","gdata"))
#####################################################
# name of the general input data file
Input_1_name= file.choose("Data") # name of the input data table as text file
# set working directory
setwd(dirname(Input_1_name))
# names of output files
# The name of the output AIC table
Output_name_1 = paste("AIC_Tab.txt")