Appendix B ! This program generates the data sets according to the simulation ! plan explained in the paper. ! The program is written in Microsoft Developer Studio (Fortran ! Power Station 4.0). ! ! Defining of variables. The variables COMMUNITY and GAMMA with ! one dimenSion is for third level, variables FAMILY and ETA with ! two dimensions are for second level, and variables CHILD, ! RESPONSE, W, and E are for the first level. ! INTEGER i,j,k, NOUT, NR,COMMUNITY(20),FAMILY(20,20) INTEGER IR,CHILD(20,20,20),RESPONSE(20,20,20) REAL RNNOF,E(20,20,20),ETA(20,20),GAMMA(20),UNIFORM,UNI REAL W(20,20,20) EXTERNAL RNNOF, UMACH,RNBIN ! ! Three files COMMUNITY, FAMILY, AND CHILD are used to save the ! sumulated data at levels 3, 2, and ! 1 respectively. The file RESPONSE is for saving the response ! at the third level. The structure of printing in the following ! files are set so that they are readable in WinBugs31. ! OPEN(1,FILE='COMMUNITY.TXT') OPEN(2,FILE='FAMILY.TXT') OPEN(3,FILE='CHILD.TXT') OPEN(4,FILE='RESPONSE.TXT') WRITE(1,9991) WRITE(2,9992) WRITE(3,9993) WRITE(4,9994) ! ! NR=20 is the sample size in each level. ! NR=20 ! ! The following three nested loops are used to generate random ! effect and a fixed effect at level 3, a random effect and a ! fixed effect at level 2, and a fixed effect and a binomial ! response at level 1. ! CALL UMACH (2, NOUT) DO 30 I=1, NR ! ! Generating Normal random number for third level as random ! effect. ! GAMMA(I)= RNNOF() ! ! Generating Binomial random number in third level for using as ! two level factor. ! CALL RNBIN (1, 1, 0.5, IR) COMMUNITY(I)=IR DO 20 j=1, NR ! ! Generating Normal random number for second level as random ! effect. ! ETA(I,J)= RNNOF() ! ! Generating Binomial random number in second level for using ! as two level factor. ! CALL RNBIN (1, 1, 0.5, IR) FAMILY(I,J)=IR DO 10 K=1, NR ! ! Generating Logistic random number for first level to establish ! logit model. ! CALL RNUN (1, UNIFORM) UNI=(UNIFORM)/(1-UNIFORM) E(I,J,K)=LOG(UNI) ! ! Generating Binomial random number in first level for using as ! two level factor. ! CALL RNBIN (1, 1, 0.5, IR) CHILD(I,J,K)=IR ! ! Generating Binomial random number in first level for using as ! response. ! W(I,J,K) is used as a latent variable to generate the binomial ! response. ! W(I,J,K)=COMMUNITY(I)+GAMMA(I)+FAMILY(I,J)+ETA(I,J) +CHILD(I,J,K)+E(I,J,K) IF(W(I,J,K).GE.0.0)THEN RESPONSE(I,J,K)=1 ELSE RESPONSE(I,J,K)=0 END IF 10 CONTINUE ! ! The following commands are provided to produce Splus format ! output readable in WinBugs31. ! For information about the splus format one may see the splus ! manual or WinBugs31 manual. ! IF(I.EQ.NR.AND.J.EQ.NR)THEN WRITE(4,9986)(RESPONSE(I,J,K),K=1,NR) WRITE(3,9986)(CHILD(I,J,K),K=1,NR) ELSE WRITE (4,9995)(RESPONSE(I,J,K),K=1,NR) WRITE (3,9996)(CHILD(I,J,K),K=1,NR) END IF 20 CONTINUE IF(I.EQ.NR)THEN WRITE(2,9986)(FAMILY(I,J),J=1,NR) ELSE WRITE (2,9998)(FAMILY(I,J),J=1,NR) END IF 30 CONTINUE WRITE (1,9997)(COMMUNITY(I),I=1,NR) WRITE(1,9990) WRITE(2,9989) WRITE(3,9988) WRITE(4,9987) 9998 FORMAT(55X,I1,19(","I1),",") 9997 FORMAT(55X,I1,19(","I1)) 9996 FORMAT(55X,I1,19(","I1),",") 9995 FORMAT(55X,I1,19(","I1),",") 9994 FORMAT("list(response=structure( .Data=c(") 9993 FORMAT("list(child=structure( .Data=c(") 9992 FORMAT("list(family=structure( .Data=c(") 9991 FORMAT("list(N=20,NUM=1,community=c(") 9990 FORMAT("))") 9989 FORMAT(".Dim=c(20,20)))") 9988 FORMAT(".Dim=c(20,20,20)))") 9987 FORMAT(".Dim=c(20,20,20)))") 9986 FORMAT(55X,I1,19(","I1),"),") END ********************************************************************** ! Appendix C ! This program generates the data sets according to the simulation ! plan explained in paper. ! Six Guassian Quadrature points are used to estimate the ! parameters in equation 9 in paper. ! The program is written in Microsoft Developer Studio (Fortran ! Power Station 4.0). ! ! Defining of variables. The variables COMMUNITY and GAMMA with ! one dimenSion is for third level, variables FAMILY and ETA with ! two dimensions are for second level, and variables ! CHILD, RESPONSE, W, and E are for the first level. ! INTEGER I,J,K, NOUT,COMMUNITY(20),FAMILY(20,20) INTEGER CHILD(20,20,20),RESPONSE(20,20,20) REAL RNNOF,E(20,20,20),ETA(20,20),GAMMA(20),UNIFORM,UNI REAL W(20,20,20) ! ! Introducing the subroutines that should be called and is ! external to this program. ! EXTERNAL RNNOF, UMACH,RNBIN ! ! Declaration of variables FOR MAXIMIZATION. NSAMPLE is the ! sample size. For the definition of the other parameters one ! should see subroutines BCONF and U4INF in the manual of Fortran ! Power Station 4.0. FCN is the subroutine for this program. ! INTEGER N,NSAMPLE PARAMETER (N=5, NSAMPLE=100) INTEGER IPARAM(7), IBTYPE REAL FVALUE, FSCALE, RPARAM(7), X(N), XGUESS(N),XLB(N) REAL XSCALE(N), XUB(N) ! ! Introducing the subroutines that should be called and is ! external to this program. ! EXTERNAL BCONF, FCN, U4INF ! DO 40 SAMPLE=1,NSAMPLE ! ! Defining the output files for simulated data and parameter ! estimation. Three files COMMUNITY, FAMILY, AND CHILD are used ! to save the sumulated data at levels 3, 2, and 1 respectively. ! The file RESPONSE is for saving the response at the third ! level.The file NONPAOUT is for saving the parameters estimate. ! OPEN(1,FILE='COMMUNITY.TXT') OPEN(2,FILE='FAMILY.TXT') OPEN(3,FILE='CHILD.TXT') OPEN(4,FILE='RESPONSE.TXT') OPEN(9,FILE='NONPAOUT.TXT') ! ! Generating the data. NR=20 is the sample size in each level. ! NR=20 ! ! The following three nested loops are used to generate a ! random effect and a fixed effect at level 3, a random effect ! and a fixed effect at level 2, and a fixed effect and a ! binomial response at level 1. ! CALL UMACH (2, NOUT) DO 30 I=1, NR ! ! Generating Normal random number for third level as ! random effect. ! GAMMA(I)= RNNOF() ! ! Generating Binomial random number in third level for using ! as two level factor. ! CALL RNBIN (1, 1, 0.5, IR) COMMUNITY(I)=IR DO 20 j=1, NR ! ! Generating Normal random number for second level as ! random effect. ! ETA(I,J)= RNNOF() ! ! Generating Binomial random number in second level for using ! as two level factor. ! CALL RNBIN (1, 1, 0.5, IR) FAMILY(I,J)=IR DO 10 K=1, NR ! ! Generating Logistic random number for first level to ! establish logit model. ! CALL RNUN (1, UNIFORM) UNI=(UNIFORM)/(1-UNIFORM) E(I,J,K)=LOG(UNI) ! ! Generating Binomial random number in first level for using ! as two level factor. ! CALL RNBIN (1, 1, 0.5, IR) CHILD(I,J,K)=IR ! ! Generating Binomial random number in first level for using ! as binary response.W(I,J,K) is used as a latent variable to ! generate the binomial response. ! W(I,J,K)=COMMUNITY(I)+GAMMA(I)+FAMILY(I,J)+ETA(I,J) +CHILD(I,J,K)+E(I,J,K) IF(W(I,J,K).GE.0.0)THEN RESPONSE(I,J,K)=1 ELSE RESPONSE(I,J,K)=0 END IF 10 CONTINUE WRITE(4,9986)(RESPONSE(I,J,K),K=1,20) WRITE(3,9986)(CHILD(I,J,K),K=1,20) 20 CONTINUE WRITE(2,9986)(FAMILY(I,J),J=1,20) 30 CONTINUE WRITE (1,9986)(COMMUNITY(I),I=1,20) 9986 FORMAT(20(1X,I1)) CLOSE (1) CLOSE (2) CLOSE (3) CLOSE (4) ! ! Introducing the parameters for calling subroutine BCONF. ! XGUESS introduces the intial values, XLB introduces the ! lower bounds,XUB introduces the upper bounds for the ! parameters, for XSCALE and FSCALE see the manual. ! DATA XGUESS/0.9E0, 0.7E0, 0.7E0, 0.5E0, 0.5E0/, XSCALE/1.0E0,1.0E0,1.0E0,1.0E0,1.0E0/, FSCALE/1.0E0/ DATA XLB/-10.0E0, -10.0E0, -10.0E0, 0.01E0,0.01E0/, XUB/10.0E0, 10.0E0, 10.0E0, 3.0E0,3.0E0/ ! ! Fitting model. ! All bounds are provided ! IBTYPE = 0 ! ! Default parameters are used ! IPARAM(1) = 0 ! ! Minimization by calling subroutine BCONF using initial ! guesses. ! FCN will be called from subroutine FCN, X is th vector ! of estimated ! parameters, N is the number of the parameters and for ! other parameters see the manual. ! CALL BCONF (FCN, N, XGUESS, IBTYPE, XLB, XUB, XSCALE, FSCALE, IPARAM, RPARAM, X, FVALUE) ! ! Printing the results ! WRITE (9,99999) X, FVALUE 99999 FORMAT (5(1X,F10.6),5X,F16.6) WRITE (NOUT,*)'SOLUTION # ',SAMPLE 40 CONTINUE END ! ! ! ! The following subroutine calculates the minus two ! log liklehood function that should be minimized by the ! main program. This gets N and X from main program and ! returns F as minus two log liklehood function. ! SUBROUTINE FCN (N, X, F) ! ! The following variables have the same definition as in ! the main program. ! INTEGER COMMUNITY(20),FAMILY(20,20),CHILD(20,20,20) INTEGER RESPONSE(20,20,20) ! ! NMASS1 and NMASS2 are the number of the mass points. ! INTEGER N,NMASS1,NMASS2,NR ! ! XX is the vector of the mass points and their ! probabilities. ! REAL X(N), F, XX(24) REAL LINEAR,EXPLINEAR,LOGISTIC,LOGISTIC1 REAL LOGISTIC2,LOGISTIC3,LOGLIK PARAMETER (NMASS1=6,NMASS2=6,NR=20) ! ! Reding data simulaled in main program. ! OPEN(5,FILE='COMMUNITY0.TXT') READ (5,9999)(COMMUNITY(I),I=1,20) OPEN(6,FILE='FAMILY0.TXT') READ (6,9999)((FAMILY(I,J),J=1,20),I=1,20) OPEN(7,FILE='CHILD0.TXT') READ (7,9999)(((CHILD(I,J,K),K=1,20),J=1,20),I=1,20) OPEN(8,FILE='RESPONSE0.TXT') READ (8,9999)(((RESPONSE(I,J,K),K=1,20),J=1,20),I=1,20) 9999 FORMAT(20(1X,I1)) CLOSE (5) CLOSE (6) CLOSE (7) CLOSE (8) ! ! Defining quadraturs and their probabilities. For each ! random effects we need 6 quadratures and their ! probabilities. ! XX(1)=3.3243 XX(2)=1.8892 XX(3)=0.6167 XX(4)=-0.6167 XX(5)=-1.8892 XX(6)=-3.3243 XX(7)=3.3243 XX(8)=1.8892 XX(9)=0.6167 XX(10)=-0.6167 XX(11)=-1.8892 XX(12)=-3.3243 XX(13)=0.0025 XX(14)=0.0886 XX(15)=0.4089 XX(16)=0.4089 XX(17)=0.0886 XX(18)=0.0025 XX(19)=0.0025 XX(20)=0.0886 XX(21)=0.4089 XX(22)=0.4089 XX(23)=0.0886 XX(24)=0.0025 ! ! Deviance evaluation using model 9 in paper. ! LOGLIK=0.0E0 DO 50 I=1,NR LOGISTIC3=0.0E0 DO 60 L=1,NMASS2 LOGISTIC2=1.0E0 DO 70 J=1,NR LOGISTIC1=0.0E0 DO 80 M=1,NMASS1 ! ! We define LOGISTIC=1.0E2 insted of one to avoid ! mathematical error then reduce 1842.068.74 from LOGLIK ! according to the model (equation 9) explained in paper. ! LOGISTIC=1.0E2 DO 90 K=1,NR LINEAR=X(1)*COMMUNITY(I)+X(2)*FAMILY(I,J)+X(3)*CHILD(I,J,K) +X(4)*XX(M)+X(5)*XX(NMASS1+L) EXPLINEAR=EXP(LINEAR) LOGISTIC=LOGISTIC*(EXP(LINEAR*RESPONSE(I,J,K))/(1+EXPLINEAR)) 90 CONTINUE LOGISTIC1=LOGISTIC1+LOGISTIC*XX(NMASS1+NMASS2+M) 80 CONTINUE LOGISTIC2=LOGISTIC2*LOGISTIC1 70 CONTINUE LOGISTIC3=LOGISTIC3+LOGISTIC2*XX(6+NMASS1+NMASS2+L) 60 CONTINUE ! ! Small positive number is added to LOGISTIC3 to avoid ! mathematical error when this ! parameter is close to zero. ! LOGLIK=LOGLIK+LOG(LOGISTIC3+0.000000000000001) 50 CONTINUE F=(-2)*(LOGLIK-1842.068074) RETURN END ********************************************************************** ! Appendix D ! This program generates the data sets according to the simulation plan explained in paper. ! Three mass points are used for Non-parametric approach to ! estimate the parameters in equation 9 introduced in paper. ! The program is written in Microsoft Developer Studio (Fortran ! Power Station 4.0). ! ! Defining of variables. The variables COMMUNITY and GAMMA with ! one dimenSion is for third level, variables FAMILY and ETA with ! two dimensions are for second level, and variables ! CHILD, RESPONSE, W, and E are for the first level. ! INTEGER i,j,k,NR,NOUT,COMMUNITY(20),FAMILY(20,20) INTEGER CHILD(20,20,20),RESPONSE(20,20,20) REAL RNNOF,E(20,20,20),ETA(20,20),GAMMA(20),UNIFORM,UNI REAL W(20,20,20) ! ! Introducing the subroutines that should be called and is external ! to this program. ! EXTERNAL RNNOF, UMACH,RNBIN ! ! Declaration of variables FOR MAXIMIZATION. NSAMPLE is the ! sample size. For the definition of the other parameters one ! should see subroutines LCONF and UMACH in the manual of Fortran ! Power Station 4.0. FCN is the subroutine for this program. ! INTEGER LDA, NCON, NEQ, NVAR, SAMPLE, NSAMPLE PARAMETER (NCON=2, NEQ=2, NVAR=17, LDA=NCON, NSAMPLE=100,NR=20) INTEGER IACT(100), MAXFCN, NACT REAL A(NCON,NVAR), ACC, ALAMDA(NVAR), B(NCON), OBJ, REAL SOL(NVAR), XGUESS(NVAR), XLB(NVAR), XUB(NVAR) ! ! Introducing the subroutines that should be called and is ! external to this program. ! EXTERNAL FCN, LCONF ! DO 40 SAMPLE=1,NSAMPLE ! ! Defining the output files for simulated data and parameter ! estimation. ! Three files COMMUNITY, FAMILY, AND CHILD are used to save ! the sumulated data at levels 3, 2, and 1 respectively. ! The file RESPONSE is for saving the response at the third ! level.The file NONPAOUT is for saving the parameters estimate. ! OPEN(1,FILE='COMMUNITY.TXT') OPEN(2,FILE='FAMILY.TXT') OPEN(3,FILE='CHILD.TXT') OPEN(4,FILE='RESPONSE.TXT') OPEN(9,FILE='NONPAOUT.TXT') ! ! Generating the data. NR=20 is the sample size in each level. ! ! ! The following three nested loops are used to generate a random ! effect and a fixed effect at level 3, a random effect and a ! fixed effect at level 2, and a fixed effect and a binomial ! response at level 1. ! CALL UMACH (2, NOUT) DO 30 I=1, NR ! ! Generating Normal random number for third level as ! random effect. ! GAMMA(I)= RNNOF() ! ! Generating Binomial random number in third level ! for using as two level factor. ! CALL RNBIN (1, 1, 0.5, IR) COMMUNITY(I)=IR DO 20 j=1, NR ! ! Generating Normal random number for second level as ! random effect. ! ETA(I,J)= RNNOF() ! ! Generating Binomial random number in second level for ! using as two level factor. ! CALL RNBIN (1, 1, 0.5, IR) FAMILY(I,J)=IR DO 10 K=1, NR ! ! Generating Logistic random number for first level to ! establish logit model. ! CALL RNUN (1, UNIFORM) UNI=(UNIFORM)/(1-UNIFORM) E(I,J,K)=LOG(UNI) ! ! Generating Binomial random number in first level for ! using as two level factor. ! CALL RNBIN (1, 1, 0.5, IR) CHILD(I,J,K)=IR ! ! Generating Binomial random number in first level for ! using as binary response. ! W(I,J,K) is used as a latent variable to generate the ! binomial response. ! W(I,J,K)=COMMUNITY(I)+GAMMA(I)+FAMILY(I,J)+ETA(I,J) +CHILD(I,J,K)+E(I,J,K) IF(W(I,J,K).GE.0.0)THEN RESPONSE(I,J,K)=1 ELSE RESPONSE(I,J,K)=0 END IF 10 CONTINUE WRITE(4,9986)(RESPONSE(I,J,K),K=1,NR) WRITE(3,9986)(CHILD(I,J,K),K=1,NR) 20 CONTINUE WRITE(2,9986)(FAMILY(I,J),J=1,NR) 30 CONTINUE WRITE (1,9986)(COMMUNITY(I),I=1,NR) 9986 FORMAT(20(1X,I1)) CLOSE (1) CLOSE (2) CLOSE (3) CLOSE (4) ! ! Fitting model. ! Minimization of -2LOGLIKELIHOOD. ! ! Sum of the probabilities for three masses should be one: ! X(12)+X(13) + X(14).EQ.1 and X(15)+X(16) + X(17).EQ.1 ! ! Defining bounds for the standard errors of random effects: ! 0 .LT. X(4) .LE. 5 and 0 .LT. X(5) .LE. 5 ! DATA A/0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0, 1.0,0.0,1.0,0.0,0.0,1.0,0.0,1.0,0.0,1.0/,B/1.0,1.0/ DATA XLB/-10,-10,-10,0.01,0.01,-100,-100,-100,-100,-100, -100,0.01,0.01,0.01,0.01,0.01,0.01/,XUB/10,10, 10,5,5,100,100,100,100,100,100,0.99,0.99, 0.99,0.99,0.99,0.99/,XGUESS/0.7,0.7,0.7,0.9,0.9, -1.0,0.0,1.0,-1.0,0.0,1.0,0.5,0.3,0.2,0.5,0.3,0.2/ DATA ACC/0.0/,MAXFCN/800/ CALL UMACH (2, NOUT) ! ! Minimization by calling subroutine LCONF using initial guesses. ! FCN will be called from subroutine FCN. For definition of the ! parameters see the manual. ! CALL LCONF (FCN, NVAR, NCON, NEQ, A, LDA, B, XLB, XUB, XGUESS, ACC, MAXFCN, SOL, OBJ, NACT, IACT, ALAMDA) ! ! Printing the parameters estimates. ! WRITE (NOUT,99998) 'Solution:',SAMPLE WRITE (9,99999) SOL,OBJ 99998 FORMAT (//, ' ', A, I4, I4) 99999 FORMAT (17(1X,F10.6),5X,F16.6) ! ! Setting the maximum function evaluation to 800. ! MAXFCN=800 40 CONTINUE END ! ! ! ! The following subroutine calculates the minus two log liklehood ! function that should be minimized by the main program. This gets ! N and X from main program and returns F as minus two ! log liklehood function. ! SUBROUTINE FCN (N, X, F) ! ! The following variables have the same definition as in ! the main program. ! INTEGER COMMUNITY(20),FAMILY(20,20),CHILD(20,20,20) INTEGER RESPONSE(20,20,20) ! ! NMASS1 and NMASS2 are the number of the mass points. ! INTEGER N,NMASS1,NMASS2,NR REAL X(N), F DOUBLE PRECISION LINEAR,EXPLINEAR,LOGISTIC,LOGISTIC1 DOUBLE PRECISION LOGISTIC2,LOGISTIC3,LOGLIK PARAMETER (NMASS1=3,NMASS2=3,NR=20) ! ! Reding simulated data from main program. ! OPEN(5,FILE='COMMUNITY.TXT') READ (5,9999)(COMMUNITY(I),I=1,NR) OPEN(6,FILE='FAMILY.TXT') READ (6,9999)((FAMILY(I,J),J=1,NR),I=1,NR) OPEN(7,FILE='CHILD.TXT') READ (7,9999)(((CHILD(I,J,K),K=1,NR),J=1,NR),I=1,NR) OPEN(8,FILE='RESPONSE.TXT') READ (8,9999)(((RESPONSE(I,J,K),K=1,NR),J=1,NR),I=1,NR) 9999 FORMAT(20(1X,I1)) CLOSE (5) CLOSE (6) CLOSE (7) CLOSE (8) ! ! Deviance evaluation using model 9 in paper. ! LOGLIK=0.0D0 DO 50 I=1,NR LOGISTIC3=0.0D0 DO 60 L=1,NMASS2 LOGISTIC2=1.0D0 DO 70 J=1,NR LOGISTIC1=0.0D0 DO 80 M=1,NMASS1 ! ! We define LOGISTIC=1.0D+7 insted of one to avoid mathematical ! error then reduce 6447.23826 from LOGLIK according to the model ! (equation 9) explained in paper. ! LOGISTIC=1.0D+7 DO 90 K=1,NR LINEAR=X(1)*COMMUNITY(I)+X(2)*FAMILY(I,J)+X(3)*CHILD(I,J,K) +X(4)*X(5+M)+X(5)*X(5+NMASS1+L) EXPLINEAR=EXP(LINEAR) LOGISTIC=LOGISTIC*(EXP(LINEAR*RESPONSE(I,J,K))/(1+EXPLINEAR)) 90 CONTINUE LOGISTIC1=LOGISTIC1+LOGISTIC*X(5+NMASS1+NMASS2+M) 80 CONTINUE LOGISTIC2=LOGISTIC2*LOGISTIC1 70 CONTINUE LOGISTIC3=LOGISTIC3+LOGISTIC2*X(8+NMASS1+NMASS2+L) 60 CONTINUE ! ! Small positive number is added to LOGISTIC3 to avoid mathematical ! error when this ! parameter is close to zero. ! LOGLIK=LOGLIK+DLOG(LOGISTIC3+0.000000000000001) 50 CONTINUE F=(-2)*(LOGLIK-6447.23826) RETURN END ***************************************************************** #Appendix E # Bugs Codes (Software WinBUGS31) for MCMC approach. Data are # simulated using codes in appendix B. model simulation { for (i in 1:N){ gamma[i]~dnorm(0.0,vargamma) for (j in 1:N){ eta[i,j]~dnorm(0.0,vareta) for (k in 1:N){ logit(mu[i,j,k])<- betacom*community[i]+betafam*family[i,j] +betachi*child[i,j,k]+gamma[i]+eta[i,j] response[i,j,k]~dbin(mu[i,j,k],NUM) llike[i,j,k]<-response[i,j,k]*log(mu[i,j,k]) +(1-response[i,j,k])*log(1-mu[i,j,k]) } } } betacom~dnorm(0.0,1.0E-6) betafam~dnorm(0.0,1.0E-6) betachi~dnorm(0.0,1.0E-6) vargamma~dgamma(1.0E-3,1.0E-3) vareta~dgamma(1.0E-3,1.0E-3) sigmaeta<-1.0/sqrt(vareta) sigmagamma<-1.0/sqrt(vargamma) llikelihood<-(-2)*sum(llike[,,]) }