lapply 不适用于 gls +piecewiseSEM 对象

如何解决lapply 不适用于 gls +piecewiseSEM 对象

我正在尝试使用 lapply 来运行模型（更具体地说是 piecewiseSEM 包中的路径模型）。这些路径模型使用单独的模型（来自 nlme 包），然后结合起来构建最终的路径模型。我不得不利用此 post 中的一些自定义函数来使模型工作。但是，现在当我尝试使用 lapply 创建的对象运行路径模型时，模型不会运行。然而，当不使用 lapply 时，它们运行得非常好。我想使用 lapply 因为我以后还想使用 Parlapply。这是一个可重现的示例：

这是我用于包 nlme 的自定义函数的代码：

library(nlme)
library(piecewiseSEM)
#### corHaversine - spatial correlation with haversine distance

# Calculates the geodesic distance between two points specified by radian latitude/longitude using Haversine formula.
# output in km
haversine <- function(x0,x1,y0,y1) {
  a <- sin( (y1 - y0)/2 )^2 + cos(y0) * cos(y1) * sin( (x1 - x0)/2 )^2
  v <- 2 * asin( min(1,sqrt(a) ) )
  6371 * v
}

# function to compute geodesic haversine distance given two-column matrix of longitude/latitude
# input is assumed in form decimal degrees if radians = F
# note fields::rdist.earth is more efficient
haversineDist <- function(xy,radians = F) {
  if (ncol(xy) > 2) stop("Input must have two columns (longitude and latitude)")
  if (radians == F) xy <- xy * pi/180
  hMat <- matrix(NA,ncol = nrow(xy),nrow = nrow(xy))
  for (i in 1:nrow(xy) ) {
    for (j in i:nrow(xy) ) {
      hMat[j,i] <- haversine(xy[i,1],xy[j,xy[i,2],2]) 
    }
  }
  as.dist(hMat)
}

## for most methods,machinery from corSpatial will work without modification
Initialize.corHaversine <- nlme:::Initialize.corSpatial
recalc.corHaversine <- nlme:::recalc.corSpatial
Variogram.corHaversine <- nlme:::Variogram.corSpatial
corFactor.corHaversine <- nlme:::corFactor.corSpatial
corMatrix.corHaversine <- nlme:::corMatrix.corSpatial
coef.corHaversine <- nlme:::coef.corSpatial
"coef<-.corHaversine" <- nlme:::"coef<-.corSpatial"

## Constructor for the corHaversine class
corHaversine <- function(value = numeric(0),form = ~ 1,mimic = "corSpher",nugget = FALSE,fixed = FALSE) {
  spClass <- "corHaversine"
  attr(value,"formula") <- form
  attr(value,"nugget") <- nugget
  attr(value,"fixed") <- fixed
  attr(value,"function") <- mimic
  class(value) <- c(spClass,"corStruct")
  value
}   # end corHaversine class
environment(corHaversine) <- asNamespace("nlme")

Dim.corHaversine <- function(object,groups,...) {
  if (missing(groups)) return(attr(object,"Dim"))
  val <- Dim.corStruct(object,groups)
  val[["start"]] <- c(0,cumsum(val[["len"]] * (val[["len"]] - 1)/2)[-val[["M"]]])
  ## will use third component of Dim list for spClass
  names(val)[3] <- "spClass"
  val[[3]] <- match(attr(object,"function"),c("corSpher","corExp","corGaus","corLin","corRatio"),0)
  val
}
environment(Dim.corHaversine) <- asNamespace("nlme")


## getCovariate method for corHaversine class
getCovariate.corHaversine <- function(object,form = formula(object),data) {
  if (is.null(covar <- attr(object,"covariate"))) {          # if object lacks covariate attribute
    if (missing(data)) {                                    # if object lacks data
      stop("need data to calculate covariate")
    }
    covForm <- getCovariateFormula(form)
    if (length(all.vars(covForm)) > 0) {                    # if covariate present
      if (attr(terms(covForm),"intercept") == 1) {       # if formula includes intercept
        covForm <- eval(parse(text = paste("~",deparse(covForm[[2]]),"-1",sep="")))    # remove intercept
      }
      # can only take covariates with correct names
      if (length(all.vars(covForm)) > 2) stop("corHaversine can only take two covariates,'lon' and 'lat'")
      if ( !all(all.vars(covForm) %in% c("lon","lat")) ) stop("covariates must be named 'lon' and 'lat'")
      covar <- as.data.frame(unclass(model.matrix(covForm,model.frame(covForm,data,drop.unused.levels = TRUE) ) ) )
      covar <- covar[,order(colnames(covar),decreasing = T)] # order as lon ... lat
    }
    else {
      covar <- NULL
    }
    
    if (!is.null(getGroupsFormula(form))) {                 # if groups in formula extract covar by groups
      grps <- getGroups(object,data = data)
      if (is.null(covar)) {
        covar <- lapply(split(grps,grps),function(x) as.vector(dist(1:length(x) ) ) )     # filler?
      } 
      else {
        giveDist <- function(el) {
          el <- as.matrix(el)
          if (nrow(el) > 1) as.vector(haversineDist(el))                       
          else numeric(0)
        }
        covar <- lapply(split(covar,giveDist )
      }
      covar <- covar[sapply(covar,length) > 0]  # no 1-obs groups
    } 
    else {                                  # if no groups in formula extract distance
      if (is.null(covar)) {
        covar <- as.vector(dist(1:nrow(data) ) )
      } 
      else {
        covar <- as.vector(haversineDist(as.matrix(covar) ) )
      }
    }
    if (any(unlist(covar) == 0)) {          # check that no distances are zero
      stop("cannot have zero distances in \"corHaversine\"")
    }
  }
  covar
}   # end method getCovariate
environment(getCovariate.corHaversine) <- asNamespace("nlme")

这是 mtcars 数据集的可重现示例/问题：

set.seed(42)  ## for sake of reproducibility
mtcars <- within(mtcars,{
  lon <- runif(nrow(mtcars))
  lat <- runif(nrow(mtcars))
})

#this makes a list of dataframes
empty_list<-replicate(n = 10,expr = mtcars,simplify = F)



#doing it the lapply method
model1<-lapply(empty_list,FUN = function(i)
  nlme::gls(disp ~  wt,correlation = corHaversine(form=~lon+lat,mimic="corSpher"),data = i)
)
model2<-lapply(empty_list,FUN = function(i)
  nlme::gls(wt ~  hp,data = i)
)

model1.2<-psem(model1[[1]],model2[[1]],data = empty_list[[1]])
summary(model1.2,.progressBar = F,standardize = "scale")

这会导致此错误：

Error in max(sapply(nm[dfdetect],nrow)) : 
  invalid 'type' (list) of argument

但是当我在没有 lapply 的情况下执行此操作时，效果很好：

model3<-nlme::gls(disp ~  wt,data = empty_list[[1]])
model4<-nlme::gls(wt ~  hp,data = empty_list[[1]])


model3.4<-psem(model3,model4)
summary(model3.4,standardize = "scale")