0

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)
Raphael
  • 1,143
  • 1
  • 11
  • 16
  • Can you add an example of your function inputs? I test it using `fm1 <- lmList(distance ~ age | Subject, Orthodont)`, but it gives me error. – agstudy Jun 21 '13 at 01:15

1 Answers1

2

I think you're looking to import dependencies into your package name space, the cheap way of doing this being to add

Imports: lme4

to your package DESCRIPTION file, and

import(lme4)

to your package NAMESPACE file. The more careful way of doing this is to enumerate all of your dependencies as described in Writing R Extensions.

Hmm you're using roxygen2; and I'm not sure how to specify imports appropriately in that context.

Martin Morgan
  • 45,935
  • 7
  • 84
  • 112
  • Disclaimer: I know very little about packaging (yet). Would simply adding `requires(lme4)` within the function work? – Scott Ritchie Jun 21 '13 at 01:33
  • Yes, but a name space is better because it guarantees that the intended functions are found (adding `require(lme4)` adds the package to the `search()` path, so someone (you) could come along and add a function that masks the intended function in lme4). – Martin Morgan Jun 21 '13 at 02:20
  • Hi Martin Morgan, Thanks so much for this very helpful information/suggestion!! I followed your advice and added lme4 both to the DESCRIPTION and NAMESPACE files and it works perfectly! In roxygen2 adding anything to the NAMESPACE is very easy: you just have to add “#’ @import lme4” on top of the specific function and it gets automatically written to the NAMESPACE. – Raphael Jun 21 '13 at 22:15
  • I also read through “Writing R Extensions.” However, since I really don’t know what specific method or class from the lme4 package is required for extracting fixed effects and standard errors from the lmer model output object, I import the full package to be on the save side, most likely it is the “mer” class but I would have to check on this. Also, while reading “Writing R Extensions” I became a little bit confused about when to list a package dependence under the “Imports” vs. the “Depends” heading in the DESCRIPTION file. Anyways, that is a minor thing since now everything works just fine. – Raphael Jun 21 '13 at 22:16