了解函数在R

如何解决了解函数在R

为了正确地适应R中的模型，我尝试了解函数的实际作用。但是，存在一些限制...

我想使用colext包中的unmarked函数来拟合回归模型。我有来自200个研究地块的鸟类物种的存在/不存在数据，这些数据分别在两年中和几个协变量中进行了调查。以下是一些示例数据：

set.seed (1)
depend1 <- sample(c(0,1),replace=TRUE,size=200)
depend2 <- sample(c(1,0),size=200)

var_g <- rnorm (200,mean = 20)*(depend + 5)
var_t1 <- var_g * (runif (200,min = 0.1,max = 0.3)+1)
var_t2 <- var_g * (runif (200,min = 0.2,max = 0.4)+1)

library (unmarked)

M <- 200 # number of Sites
Times <- 2 # num primary sample periods

occurences <- array (data = c (depend1,depend2),dim = c (M,Times)) %>% as.matrix()
covariatesSite <- data.frame (var_g = var_g)
covariatesSiteYear <- list (var_t = data.frame (year1 = var_t1,year2 = var_t2))

my_data <- unmarkedMultFrame (y = occurences,numPrimary = Times,siteCovs = covariatesSite,yearlySiteCovs = covariatesSiteYear
                                )

我尝试通过colext函数通过动态占用模型来模拟物种的存在/不存在。空模型看起来像这样：

m0 <- colext (psiformula= ~  1,gammaformula = ~  1,epsilonformula = ~ 1,pformula = ~ 1,data = my_data,method="BFGS")
summary(m0)

使用psi公式，可以对“占用率”进行建模，例如：

m1 <- colext (psiformula= ~  var_g  + 1,method="BFGS")
summary(m1)

但是，有

一个警告：对于这种模型类型，这四个参数之一必须是 派生也就是说，您只能获取 研究的第一个季节及其后的模式必须推断出来。

我在https://darinjmcneil.weebly.com/multi-season-occupancy.html处找到了这句话，上面也写着：

因为在多季节占用模型中需要关注的参数 通常是γ和ε，通常是psi 参数。这意味着当我们估计“ psi”或 “入住率”应被视为“第一季的入住率”。 必须从殖民化和 灭绝。

现在让我们假设var_g随时间推移不稳定。但是，我只有时间实例2上的var_g数据。因此，我认为“扭转”一切都是一个好主意：第二年适合psi，第一年不适合。我的方法是更改​​my_data的结构：

occurences <- array (data = c (depend2,depend1),Times)) %>% as.matrix()
covariatesSite <- data.frame (var_g = var_g)
covariatesSiteYear <- list (var_t = data.frame (year2 = var_t2,year1 = var_t1))

my_data <- unmarkedMultFrame (y = occurences,yearlySiteCovs = covariatesSiteYear
                                )

，然后使用与以前相同的模型：

m1 <- colext (psiformula= ~  var_g  + 1,method="BFGS")
summary(m1)

现在，我认为伽玛不再灭绝，而是定居（反之亦然）。而且我认为psi的估算是指第二年。但是，我对此不太确定。所以我想知道函数的作用。我尝试过

> colext
function (psiformula = ~1,gammaformula = ~1,epsilonformula = ~1,pformula = ~1,data,starts,method = "BFGS",se = TRUE,...) 
{
    K <- 1
    data@y[data@y > K] <- K
    y <- getY(data)
    J <- numY(data)/data@numPrimary
    M <- nrow(y)
    nY <- ncol(y)/J
    n.det <- sum(apply(y > 0,1,any,na.rm = TRUE))
    fc <- match.call()
    fc[[1]] <- as.name("colext.fit")
    formula <- list(psiformula = psiformula,gammaformula = gammaformula,epsilonformula = epsilonformula,pformula = pformula)
    fc$formula <- as.name("formula")
    fc$bootstrap.se <- fc$covdata.site <- fc$covdata.obs <- fc$data <- fc$B <- fc$psiformula <- fc$gammaformula <- fc$epsilonformula <- fc$pformula <- NULL
    fc$data <- as.name("data")
    fc$J <- as.name("J")
    fc$method <- as.name("method")
    fc$getHessian <- as.name("se")
    fc$se <- NULL
    if (missing(starts)) {
        fc$starts <- NULL
    }
    else {
        fc$starts <- eval(starts)
    }
    extras <- match.call(expand.dots = FALSE)$...
    if (length(extras) > 0) {
        existing <- !is.na(match(names(extras),names(fc)))
        for (a in names(extras)[existing]) fc[[a]] <- extras[[a]]
        if (any(!existing)) {
            fc <- as.call(c(as.list(fc),extras[!existing]))
        }
    }
    fm <- eval(fc)
    fm$n.det <- n.det
    opt <- fm$opt
    nP <- fm$nP
    M <- fm$M
    nDP <- fm$nDP
    nGP <- fm$nGP
    nEP <- fm$nEP
    nSP <- fm$nSP
    covMat <- invertHessian(opt,nP,se)
    ests <- opt$par
    names(ests) <- fm$mle$names
    fmAIC <- 2 * opt$value + 2 * nP
    psiParams <- ests[1:nSP]
    colParams <- ests[(nSP + 1):(nSP + nGP)]
    extParams <- ests[(nSP + nGP + 1):(nSP + nGP + nEP)]
    detParams <- ests[(nSP + nGP + nEP + 1):nP]
    psi <- unmarkedEstimate(name = "Initial",short.name = "psi",estimates = psiParams,covMat = as.matrix(covMat[1:nSP,1:nSP]),invlink = "logistic",invlinkGrad = "logistic.grad")
    col <- unmarkedEstimate(name = "Colonization",short.name = "col",estimates = colParams,covMat = as.matrix(covMat[(nSP + 
            1):(nSP + nGP),(nSP + 1):(nSP + nGP)]),invlinkGrad = "logistic.grad")
    ext <- unmarkedEstimate(name = "Extinction",short.name = "ext",estimates = extParams,covMat = as.matrix(covMat[(nSP + 
            nGP + 1):(nSP + nGP + nEP),(nSP + nGP + 1):(nSP + 
            nGP + nEP)]),invlinkGrad = "logistic.grad")
    det <- unmarkedEstimate(name = "Detection",short.name = "p",estimates = detParams,covMat = as.matrix(covMat[(nSP + 
            nGP + nEP + 1):nP,(nSP + nGP + nEP + 1):nP]),invlinkGrad = "logistic.grad")
    estimateList <- unmarkedEstimateList(list(psi = psi,col = col,ext = ext,det = det))
    umfit <- new("unmarkedFitColExt",fitType = "colext",call = match.call(),formula = as.formula(paste(unlist(formula),collapse = " ")),psiformula = psiformula,gamformula = gammaformula,epsformula = epsilonformula,detformula = pformula,data = data,sitesRemoved = fm$designMats$removed.sites,estimates = estimateList,AIC = fmAIC,opt = opt,negLogLike = opt$value,nllFun = fm$nll,projected = fm$projected,projected.mean = fm$projected.mean,smoothed = fm$smoothed,smoothed.mean = fm$smoothed.mean)
    return(umfit)
}
<bytecode: 0x0000000020c16f38>
<environment: namespace:unmarked>

但是，我无法确定psi的估算方式和估算位置（尤其是现在估算是否涉及第二年）。有办法以某种方式更深入地研究该功能吗？