- Notifications
You must be signed in to change notification settings - Fork9
Open
Description
Hi Ken,
I'm not sure whether the use of a bayes R2 would be broadly suitable for occupancy models given the limitations with its application in logistic regression and often low values. However I had a crack at getting a bayes R2 for the state, detection submodels (and their combined value) using the residual method (https://avehtari.github.io/bayes_R2/bayes_R2.html). A bit of a messy function and example below.
# R2 example for logistic occupancy model library(ubms)#> Loading required package: unmarked#> Loading required package: lattice#>#> Attaching package: 'ubms'#> The following objects are masked from 'package:unmarked':#>#> fitList, projected#> The following object is masked from 'package:base':#>#> gamma# R2 functionbayes_R2_res_ubms<-function(fit,re.form=NULL,draws=draws) {# Get the observed occupancy at each site for each observation periody<-ubms::getY(fit)# Get the observed occupancy at each sitey_mod<-matrix(apply(y,1,FUN=function(x) max(x,na.rm=TRUE),simplify=TRUE),ncol=1)# Get the linear predictor for the 'state'ypred<-ubms::posterior_linpred(fit,transform=TRUE,submodel="state",re.form=re.form,draws=draws)yp<-ubms::posterior_linpred(fit,transform=TRUE,submodel="det",re.form=re.form,draws=draws)# Assume perfect detection (no effect of detection submodel)yp1<-ypyp1[!is.na(yp1)]<-1# For each draw obtain the probability for each site that a detection will occur (0-1).# Probably a less messy way to do thisys<-list("state"=yp,"detection"=yp1)ypred_mod<-list()r2<-list()for(jin1:length(ys)) {# detection * stateypred_prob<-ys[[j]]# draws x sites x obsypred_cond<-array(data=NA,dim= c(draws, dim(ypred)[2],fit@response@max_obs))ypred_mod[[j]]<-ypred# loop over drawsfor(iin1:draws) {# repeat the latent state by observation periodsypred_vec<- rep(ypred[i,],each=fit@response@max_obs)# det * stateypred_prob[i,]<-ypred_prob[i,]*ypred_vec# force into matrix with nsites X nobsypred_cond[i,,]<-matrix(ypred_prob[i,],ncol=fit@response@max_obs,byrow=TRUE)# 1 minus the probability that all observations are zero = at least 1 detectionypred_mod[[j]][i,]<-1-apply(1-ypred_cond[i,,],1,function(x) { prod(x,na.rm=TRUE) }) }if (fit@response@z_dist=="binomial"&& NCOL(y_mod)==2) {trials<- rowSums(y_mod)y_mod<-y_mod[,1]ypred_mod[[j]]<-ypred_mod[[j]]%*% diag(trials) }e<--1* sweep(ypred_mod[[j]],2,y_mod)var_ypred<- apply(ypred_mod[[j]],1,var)var_e<- apply(e,1,var)r2[[j]]<-var_ypred/ (var_ypred+var_e) }r2[[3]]<-r2[[1]]-r2[[2]] names(r2)<- c("both","state","det")return(r2) }# Data simulation set.seed(123)dat_occ<-data.frame(x1=rnorm(500))dat_p<-data.frame(x2=rnorm(500*5))y<-matrix(NA,500,5)z<- rep(NA,500)b<- c(0.4,-0.5,0.3,0.5)re_fac<-factor(sample(letters[1:26],500,replace=T))dat_occ$group<-re_facre<- rnorm(26,0,1.2)re_idx<- as.numeric(re_fac)idx<-1for (iin1:500){z[i]<- rbinom(1,1, plogis(b[1]+b[2]*dat_occ$x1[i]+re[re_idx[i]]))for (jin1:5){y[i,j]<-z[i]*rbinom(1,1, plogis(b[3]+b[4]*dat_p$x2[idx]))idx<-idx+1 } }# unmarked frameumf<- unmarkedFrameOccu(y=y,siteCovs=dat_occ,obsCovs=dat_p)# model options(mc.cores=3)#number of parallel cores to usefm<- stan_occu(~x2~x1+ (1|group),umf,chains=3) bayes_R2_res_ubms(fm,draws=100)#> $both#> [1] 0.2436708 0.1813600 0.1861829 0.1968885 0.2093023 0.2065249 0.2262449#> [8] 0.1605566 0.2356042 0.2089506 0.2313770 0.2569806 0.2476316 0.2080628#> [15] 0.1884270 0.2458283 0.2353730 0.2410739 0.2087004 0.2352926 0.2182499#> [22] 0.2195304 0.2197247 0.2299345 0.1642092 0.2512100 0.2537433 0.1866575#> [29] 0.2420921 0.1616609 0.1850477 0.2040931 0.2660747 0.1736265 0.2113426#> [36] 0.2345318 0.2077254 0.1940378 0.2166815 0.2828791 0.2590453 0.1380645#> [43] 0.2494774 0.2094812 0.1596296 0.2125061 0.2406926 0.2306142 0.2238410#> [50] 0.1424783 0.2422054 0.2399446 0.2433104 0.2083924 0.2163885 0.2143759#> [57] 0.2868863 0.2545129 0.2662583 0.2061908 0.2272771 0.1880253 0.2036795#> [64] 0.2259326 0.2396712 0.1560981 0.2029403 0.2575136 0.2447248 0.2600542#> [71] 0.1503401 0.1612022 0.2184497 0.2333725 0.2235703 0.2745752 0.1902804#> [78] 0.2021373 0.1704586 0.2658361 0.2470988 0.2298406 0.2281244 0.2116103#> [85] 0.2407505 0.2564699 0.2642025 0.2530184 0.2964650 0.1884810 0.1804507#> [92] 0.2371928 0.2304418 0.2121007 0.2515888 0.2200955 0.1956262 0.2275631#> [99] 0.2210517 0.2139585#>#> $state#> [1] 0.17605019 0.10263418 0.10580724 0.13545364 0.14482887 0.13800919#> [7] 0.17355460 0.08628623 0.17066872 0.14583832 0.16249959 0.20375155#> [13] 0.18365678 0.13735974 0.11535767 0.17160635 0.16365889 0.18999839#> [19] 0.16342425 0.17257228 0.16156030 0.15115874 0.15180402 0.17367579#> [25] 0.09177282 0.18924085 0.20799994 0.13745090 0.18211136 0.09588723#> [31] 0.11718993 0.14892238 0.22422619 0.13182738 0.15381415 0.17530606#> [37] 0.14920495 0.12388629 0.15416240 0.24336614 0.20041741 0.08053604#> [43] 0.20827089 0.14264055 0.09731327 0.14357953 0.18287715 0.16933439#> [49] 0.15418602 0.07828137 0.17981013 0.17819720 0.19682582 0.15552595#> [55] 0.15226781 0.15547727 0.24004102 0.20505126 0.22142930 0.14618668#> [61] 0.15626259 0.12546194 0.14349652 0.16704033 0.18849389 0.08944449#> [67] 0.12046083 0.22172548 0.19268544 0.20941069 0.07355600 0.11754425#> [73] 0.15150645 0.17893636 0.16290501 0.22485436 0.12711429 0.15046677#> [79] 0.11693493 0.22099586 0.18770485 0.17481347 0.16442692 0.14398650#> [85] 0.17886823 0.20548550 0.20694618 0.20829883 0.24947301 0.11531789#> [91] 0.10101353 0.18807278 0.18736563 0.14776186 0.19111258 0.16002263#> [97] 0.11816607 0.17083967 0.16663188 0.15840369#>#> $det#> [1] 0.06762062 0.07872586 0.08037562 0.06143486 0.06447340 0.06851572#> [7] 0.05269027 0.07427036 0.06493549 0.06311229 0.06887738 0.05322905#> [13] 0.06397480 0.07070302 0.07306933 0.07422198 0.07171414 0.05107554#> [19] 0.04527614 0.06272035 0.05668959 0.06837161 0.06792069 0.05625875#> [25] 0.07243638 0.06196911 0.04574332 0.04920663 0.05998069 0.06577362#> [31] 0.06785781 0.05517070 0.04184849 0.04179915 0.05752846 0.05922575#> [37] 0.05852047 0.07015146 0.06251909 0.03951300 0.05862791 0.05752848#> [43] 0.04120654 0.06684066 0.06231634 0.06892656 0.05781542 0.06127982#> [49] 0.06965496 0.06419698 0.06239527 0.06174742 0.04648460 0.05286642#> [55] 0.06412072 0.05889858 0.04684531 0.04946166 0.04482896 0.06000413#> [61] 0.07101450 0.06256336 0.06018296 0.05889223 0.05117732 0.06665366#> [67] 0.08247942 0.03578808 0.05203937 0.05064353 0.07678410 0.04365799#> [73] 0.06694323 0.05443615 0.06066532 0.04972085 0.06316610 0.05167058#> [79] 0.05352365 0.04484027 0.05939398 0.05502711 0.06369745 0.06762376#> [85] 0.06188222 0.05098439 0.05725632 0.04471958 0.04699198 0.07316309#> [91] 0.07943715 0.04912006 0.04307613 0.06433879 0.06047627 0.06007289#> [97] 0.07746012 0.05672347 0.05441978 0.05555479
Created on 2021-11-11 by thereprex package (v2.0.1)
Metadata
Metadata
Assignees
Labels
No labels