CDOoDocuments.StdDocumentDescDocuments.DocumentDescContainers.ViewDescViews.ViewDescStores.StoreDescDocuments.ModelDescContainers.ModelDescModels.ModelDescStores.ElemDesc TextViews.StdViewDescTextViews.ViewDescTextModels.StdModelDescTextModels.ModelDesc TextModels.AttributesDesc1$Courier New  1$Courier New CourierPS -y  #------------------------------------------------------------------------------------------------------------------------------ Version 1 #------------------------------------------------------------------------------------------------------------------------------ model{ # -------------------------------------------------------------------- # model's likelihood for (i in 1:nf){ # scored goals for ( j in 1:fx[i] ){ z[i,j] <- y[i] z[i,j]~dpois( lambda ) } } # # log-normal prior lambda ~ dlnorm( 0.0, 0.001) for (i in 1:n.rep){y.rep[i]~dpois( lambda ) } for (k in 1:nf){ for (i in 1:n.rep){ freq.bin[k,i] <- equals( y.rep[i], y[k] ) } freq[k] <- sum( freq.bin[k,1:n.rep] ) rel.freq[k] <- freq[k]/n.rep } # # -------------------------------------------------------------------- # calculation of ppo's for (i in 1:nf){ # posterior predictive ordinate ppo[i] <- exp( -lambda + y[i]*log(lambda) - logfact(y[i]) ) # inverse of the conditional predictive ordinate icpo[i] <- 1/ppo[i] } # -------------------------------------------------------------------- # log-score ls <- - inprod( fx[], ppo[] ) # -------------------------------------------------------------------- # chisq for (i in 1:n.rep){chisq.rep.vec[i] <- pow(y.rep[i]-lambda,2)/lambda} chisq.rep <- sum( chisq.rep.vec[] ) for (k in 1:nf){ chisq.obs.vec[k] <- pow(y[k]-lambda,2)/lambda } chisq.obs <- inprod( fx[], chisq.obs.vec[] ) chisq.diff <- chisq.rep - chisq.obs chisq.p <- step(chisq.diff) # # -------------------------------------------------------------------- # deviance for (i in 1:n.rep){ dev.rep.vec[i] <- -lambda + y.rep[i]*log(lambda) - logfact(y.rep[i]) } dev.rep <- -2*sum( dev.rep.vec[] ) for (k in 1:nf){ dev.obs.vec[k] <- -lambda + y[k]*log(lambda) - logfact(y[k]) } dev.obs <- -2*inprod( fx[], dev.obs.vec[] ) dev.diff <- dev.rep - dev.obs dev.p <- step(dev.diff) # ------------------------------------------------- # chisq based on frequencies for(i in 1:nf){ f[i] <- exp( -lambda + y[i]*log(lambda) - logfact(y[i]) ) chisq.freq.obs.vec[i]<- pow(fx[i]/n - f[i],2)/(f[i]*(1-f[i])) chisq.freq.rep.vec[i]<- pow(rel.freq[i] - f[i],2)/(f[i]*(1-f[i])) } chisq.freq.obs<- n*sum(chisq.freq.obs.vec[]) chisq.freq.rep<- n*sum(chisq.freq.rep.vec[]) chisq.freq.p <-step( chisq.freq.rep-chisq.freq.obs ) # -------------------------------------------------------------------- # checking over/under dispersion # calculation of replicated DI di.rep <- mean( y.rep[] )/ pow( sd(y.rep[]), 2) # calculation of observed di (tabulated data) n <- sum(fx[]) # sample size bar.y <- inprod( fx[], y[] )/n # sample mean # sample variance components for (i in 1:nf){ ss[i]<- fx[i] * pow( y[i] - bar.y, 2 ) } var.y <- sum(ss[])/(n-1) # sample variance di.obs <- bar.y/var.y # observed DI di.diff <- di.rep - di.obs # DI difference di.p <- step(di.diff) # DI p-value } INITS list( lambda=1.0 ) DATA (LIST) list( nf=6, n.rep=19,y=c(0,1,2,3,4,5), fx=c(2, 3,4,6,3,1) ) # scored list( nf=6, n.rep=19, y=c(0,1,2,3,4,5), fx=c(8,10,1,0,0,0) ) # conceded TextControllers.StdCtrlDescTextControllers.ControllerDescContainers.ControllerDescControllers.ControllerDesc TextRulers.StdRulerDescTextRulers.RulerDescTextRulers.StdStyleDescTextRulers.StyleDescZTextRulers.AttributesDesc$ ZGo * ,[ @Documents.ControllerDesc t]s ' `h*