I have lately tried to organize all my functions in packages to avoid having hundreds of lines of code on top of each script. I have written a number of functions that I use after running a set of LME4 multilevel models. These functions are designed to generate nice output tables (I know there are packages available that do this, but these are usually not flexible enough to customize the tables to my needs). Here is an example of one function that takes a list of lmer models (e.g., fm0, fm1, fm2) and combines the fixed effects parameters and their respective standard errors in one output table (which I use later on together with other model statistics to generate a customized regression table).
#' @rdname f.regTableFix1
#' @title Fixed effects table for lmer models
#' @description This function produces a regression table for multiple models of type lmer (lme4).
#' @param mList list of models to be used for the table
#' @return The output contains fixed effects (b,Std.Error) parameter.
#' @export
#'
f.regTableFix1=function(mList){
# Obtain fixed effects parameter
for(m in 1:length(mList)){
#m=2
# Declare variables
fix=fixef(mList[[m]]) #obtains fixed effect estimates
stder=sqrt(diag(vcov(mList[[m]]))) #obtains respective standard errors
if(m==1){
masterFix=data.frame(variables=as.character(names(fix)),b=fix,se=stder,row.names=NULL,stringsAsFactors=F)
names(masterFix)[!names(masterFix)=="variables"]=paste(names(masterFix)[!names(masterFix)=="variables"],m,sep=".")
}else{
addFix=data.frame(variables=as.character(names(fix)),b=fix,se=stder,row.names=NULL,stringsAsFactors=F)
names(addFix)[!names(addFix)=="variables"]=paste(names(addFix)[!names(addFix)=="variables"],m,sep=".")
masterFix=merge(masterFix,addFix,by="variables",all=T,sort=F)
}
}
return(masterFix)
}
If I just use the function in my script (defined on top), it works fine. However, after I include the function in a package, it throws the following error (I generate the package in the old fashion way by using a system call to R CMD build and install the package using R CMD INSTALL).
> tableFix=f.regTableFix1(modList)
Error in as.integer(x) :
cannot coerce type 'S4' to vector of type 'integer'
This error seems to be somehow connected to the use of an S4 object in my function. Unfortunately, the fixed effects parameters and standard errors come from the lmer model (S4 object) and I can’t change this. I have tried multiple different ways of obtaining the fixed effect parameters and standard errors including the following:
coef(summary(mList[[m]]))
summary(mList[[m]])@coefs
Using these different specifications all work fine if I use them in a function, declared at the top of my script, but none of them works if I put the function into a package. Only the error message changes. For example, I get the error message “Error: $ operator is invalid for atomic vectors” when using “coef(summary(mList[[m]]))”.
All other functions in my package work fine, so I guess it is not a general issue with the way of how I create the package. Has anyone experienced this type of problems before? Are there any suggestions of how to successfully use S4 objects as input for functions in packages? Thanks so much for any help/suggestions/comments!
Best,
Raphael
Edits:
Sorry for not providing some example data in the first place. But here it is now:
library(lme4)
# Prepare some example data
edata=Dyestuff #data automatically available through lme4
edata$var1=rnorm(30)
edata$var2=rnorm(30)
edata$var3=rnorm(30)
# Run multilevel models
fm0=lmer(Yield~1+(1|Batch), data=edata,
na.action=na.omit, family="gaussian", verbose=F)
fm1=lmer(Yield~1+var1+(1|Batch), data=edata,
na.action=na.omit, family="gaussian", verbose=F)
fm2=lmer(Yield~1+var1+var2+var3+(1|Batch), data=edata,
na.action=na.omit, family="gaussian", verbose=F)
# Store model outputs in list
modList=list(fm0, fm1, fm2)
# Obtain fixed effects output table with above discussed function
fixTable=f.regTableFix1(modList)