CDOoDocuments.StdDocumentDescDocuments.DocumentDescContainers.ViewDescViews.ViewDescStores.StoreDesc Documents.ModelDescContainers.ModelDescModels.ModelDescStores.ElemDesc* " TextViews.StdViewDescTextViews.ViewDesc TextModels.StdModelDescTextModels.ModelDesc  TextModels.AttributesDesc'* * Az# # chapter 9 - example 6 # logistic regression model for 2x2 table with correlated binary responses # model 3 - fixed p for each repeated measurement model{ pi <- 3.14 for (k in 1:4){ for (i in (F[k]+1):F[k+1]){ # ------------------------------------------------------------------------- # model for the 1st individual (Treatment A) of pair i # ------------------------------------------------------------------------- y1[i]<-Y1[ k ] # definition of y1 y1[i]~dbern( p1[i] ) # stochastic component for y1 logit(p1[i])<- beta[1] # linear predictor and link for y1 # ------------------------------------------------------------------------- # model for the 2nd individual (Treatment B) of pair i # ------------------------------------------------------------------------- y2[i]<-Y2[ k] # definition of y2 y2[i]~dbern( p2[i] ) # stochastic component for y2 logit(p2[i])<- beta[1]+beta[2] # linear predictor and link for y2 # ------------------------------------------------------------------------- # random effects for individuals belonging in i pair # ------------------------------------------------------------------------- # # squared residuals sq.res1[i] <- pow(y1[i]-p1[i],2) sq.res2[i] <- pow(y2[i]-p2[i],2) }} # ------------------------------------------------------------------------- # PRIORS # ------------------------------------------------------------------------- beta[1]~dnorm(0,0.001) # prior for beta[1] (fixed effect) beta[2]~dnorm(0,0.001) # prior for beta[1] (fixed effect) # ------------------------------------------------------------------------- # LOGICAL NODES # ------------------------------------------------------------------------- # # estimating variance components # ss.res <- sum(sq.res1[1:n]) + sum(sq.res2[1:n]) overallp <- (sum(y1[1:n]) + sum(y2[1:n]))/(2*n) for (i in 1:n){ ssy1[i] <- pow( y1[i]-overallp, 2) ssy2[i] <- pow( y2[i]-overallp, 2) } ss.y <- sum( ssy1[1:n] ) + sum( ssy2[1:n] ) R2bin <- 1 - ss.res/ss.y } INITS list( beta=c(0,0)) DATA list( n=60, F=c(0,7, 9,23,60), Y1=c(0,1,0,1), Y2=c(0,0,1,1) ) # # F are accumulated counts # actual frequency data are f=c(7,2,14,37) AGRESTI'S DATA - RATINGS OF PRIME MINISTER - Table 12.1, page 494 list( n=1600, F=c(0, 794, 880, 1030, 1600), Y1=c(1,0,1,0), Y2=c(1,1,0,0) ) # actual frequency data are f=c(794,86,150,570) TextControllers.StdCtrlDescTextControllers.ControllerDescContainers.ControllerDescControllers.ControllerDesc TextRulers.StdRulerDescTextRulers.RulerDescTextRulers.StdStyleDescTextRulers.StyleDescZTextRulers.AttributesDesc$ Zo * ,[ @Documents.ControllerDesc Ws,! [h$