Add fitting function for community occupancy model (#31)

* Initial working occuComm function

* Support for some species covs

* Fix broken predict and obsNum methods

* Add unmarkedFrameOccuComm methods

* getP, fitted, residuals, SSE, getY, replaceY, simulate methods

* Plot method for unmarkedFitOccuComm

* Handle newdata in predict

* Length-S species covs

* addMean argument to randomTerms

* Fix some bugs with species covs

* Better summary method

* Fix some bugs

* Allow providing y as an array

* Don't automatically add species to predict when using newdata

* Add vignette

* Add docs

* Add tests

* Add reference to docs

* Bump version

Fixes #31

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: unmarked
-Version: 1.4.3.9008
-Date: 2025-01-28
+Version: 1.4.3.9009
+Date: 2025-02-05
 Type: Package
 Title: Models for Data from Unmarked Animals
 Authors@R: c(
@@ -63,6 +63,7 @@ Collate: 'classes.R' 'unmarkedEstimate.R' 'unmarkedFrame.R'
     'IDS.R'
     'goccu.R'
     'occuCOP.R'
+    'occuComm.R'
     'deprecated_sim_power.R'
     'RcppExports.R'
     'zzz.R'

NAMESPACE

@@ -23,7 +23,7 @@ importFrom(Rcpp, evalCpp)
 export(occu, occuFP, occuRN, pcount, pcountOpen, multinomPois, distsamp,
        colext, gmultmix, gdistsamp, gpcount, occuPEN, occuPEN_CV, occuMulti,
        occuMS, computeMPLElambda, pcount.spHDS, occuTTD, distsampOpen,
-       multmixOpen, nmixTTD, gdistremoval, goccu, occuCOP, IDS)
+       multmixOpen, nmixTTD, gdistremoval, goccu, occuCOP, IDS, occuComm)
 
 export(removalPiFun, doublePiFun)
 export(makeRemPiFun, makeCrPiFun, makeCrPiFunMb, makeCrPiFunMh)
@@ -41,7 +41,7 @@ exportMethods(backTransform, coef, confint, fitted, getData,
               "yearlySiteCovs<-", "[", smoothed, projected, nonparboot, logLik,
               ranef, bup, crossVal, posteriorSamples, sigma, randomTerms, 
               optimizePenalty, unmarkedPowerList, plotEffectsData, plotEffects,
-              getL)
+              getL, richness)
 
 S3method("print", "unmarkedPostSamples")
 
@@ -52,7 +52,7 @@ export(fitList,
     unmarkedFramePCO, unmarkedFrameGDS, unmarkedFrameGPC, unmarkedFrameOccuMulti,
     unmarkedFrameOccuMS, unmarkedFrameOccuTTD, unmarkedFrameDSO, 
     unmarkedFrameMMO, unmarkedFrameGDR, unmarkedFrameGOccu,
-    unmarkedFrameOccuCOP)
+    unmarkedFrameOccuCOP, unmarkedFrameOccuComm)
 
 # Formatting
 export(csvToUMF, formatLong, formatWide, formatMult, formatDistData)

R/mixedModelTools.R

@@ -278,7 +278,8 @@ setMethod("sigma", "unmarkedFit", function(object, type, level=0.95, ...){
 setGeneric("randomTerms", function(object, ...) standardGeneric("randomTerms"))
 
 
-setMethod("randomTerms", "unmarkedEstimate", function(object, level=0.95, ...){
+setMethod("randomTerms", "unmarkedEstimate", 
+          function(object, level=0.95, addMean = FALSE, ...){
 
   rv <- object@randomVarInfo
   if(length(rv)==0){
@@ -303,6 +304,25 @@ setMethod("randomTerms", "unmarkedEstimate", function(object, level=0.95, ...){
   b_var <- object@estimates[rv_idx]
   b_se <- sqrt(diag(object@covMat[rv_idx,rv_idx,drop=FALSE]))
 
+  if(addMean){
+    #b_se <- rep(NA, length(b_var)) 
+    fixed <- object@estimates[!rv_idx]
+
+    for (i in unique(Name)){
+      if(is.na(fixed[i])) next
+      b_var[Name == i] <- b_var[Name == i] + fixed[i] 
+
+      # Calculate SE of sum
+      fixed_ind <- which(names(fixed) == i)
+      rand_inds <- which(Name == i) + length(fixed)
+      var_fixed <- object@covMat[fixed_ind, fixed_ind]
+      vars_rand <- diag(object@covMat[rand_inds, rand_inds,drop=FALSE])
+      covvar <- object@covMat[fixed_ind, rand_inds]
+      b_se[Name == i] <- sqrt(var_fixed + vars_rand + 2 * covvar)      
+    }
+
+  }
+
   z <- qnorm((1-level)/2, lower.tail = FALSE)
   lower <- b_var - z*b_se
   upper <- b_var + z*b_se
@@ -312,10 +332,10 @@ setMethod("randomTerms", "unmarkedEstimate", function(object, level=0.95, ...){
 })
 
 
-setMethod("randomTerms", "unmarkedFit", function(object, type, level=0.95, ...){
+setMethod("randomTerms", "unmarkedFit", function(object, type, level=0.95, addMean = FALSE, ...){
 
   if(!missing(type)){
-    return(randomTerms(object[type], level))
+    return(randomTerms(object[type], level, addMean))
   }
 
   has_random <- sapply(object@estimates@estimates,
@@ -324,7 +344,7 @@ setMethod("randomTerms", "unmarkedFit", function(object, type, level=0.95, ...){
     stop("No random effects in this model", call.=FALSE)
   }
   keep <- object@estimates@estimates[has_random]
-  out <- lapply(keep, randomTerms, level=level)
+  out <- lapply(keep, randomTerms, level=level, addMean=addMean)
   out <- do.call(rbind, out)
   rownames(out) <- NULL
   out