В последнее время я пытался организовать все свои функции в пакетах, чтобы не создавать сотни строк кода поверх каждого скрипта. Я написал ряд функций, которые использую после запуска набора многоуровневых моделей LME4. Эти функции предназначены для создания хороших выходных таблиц (я знаю, что для этого есть пакеты, но они обычно недостаточно гибкие, чтобы настраивать таблицы в соответствии с моими потребностями). Вот пример одной функции, которая принимает список моделей lmer (например, fm0, fm1, fm2) и объединяет параметры фиксированных эффектов и их соответствующие стандартные ошибки в одной выходной таблице (которую я использую позже вместе со статистикой других моделей для создать настраиваемую таблицу регрессии).
#' @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)
}
Если я просто использую функцию в моем скрипте (определенном сверху), она работает нормально. Однако после включения функции в пакет возникает следующая ошибка (я генерирую пакет по старинке, используя системный вызов для сборки R CMD и устанавливаю пакет с помощью R CMD INSTALL).
> tableFix=f.regTableFix1(modList)
Error in as.integer(x) :
cannot coerce type 'S4' to vector of type 'integer'
Эта ошибка, похоже, как-то связана с использованием объекта S4 в моей функции. К сожалению, параметры фиксированных эффектов и стандартные ошибки исходят из модели lmer (объект S4), и я не могу это изменить. Я пробовал несколько разных способов получения параметров фиксированного эффекта и стандартных ошибок, включая следующие:
coef(summary(mList[[m]]))
summary(mList[[m]])@coefs
Использование этих различных спецификаций отлично работает, если я использую их в функции, объявленной в верхней части моего скрипта, но ни одна из них не работает, если я помещаю функцию в пакет. Меняется только сообщение об ошибке. Например, я получаю сообщение об ошибке «Ошибка: оператор $ недопустим для атомарных векторов» при использовании «coef (summary (mList [[m]]))». Все остальные функции в моем пакете работают нормально, поэтому я думаю, это не общая проблема, связанная с тем, как я создаю пакет. Кто-нибудь раньше сталкивался с подобными проблемами? Есть ли какие-либо предложения о том, как успешно использовать объекты S4 в качестве входных данных для функций в пакетах? Большое спасибо за любую помощь / предложения / комментарии!
Бест, Рафаэль
Изменения. Приносим извинения за то, что изначально не предоставили некоторые примеры данных. Но вот оно сейчас:
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)
fm1 <- lmList(distance ~ age | Subject, Orthodont)
, но это дает мне ошибку. - person agstudy   schedule 21.06.2013