Movatterモバイル変換


[0]ホーム

URL:


Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

Sign up
Appearance settings

bayes R2 for occupancy models #58

Open
@JustinCally

Description

@JustinCally

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

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions


      [8]ページ先頭

      ©2009-2025 Movatter.jp