/*Likelihood for the migration model Feb 19'07 added the writing of individual deviances to an external file. May 8'07 double-to-single tags was programmed incorrectly and this is corrected August 07 Data changed slightly when strata borders were double checked. 4 fish removed from dataset. Oct 07 New design matrices for the design_p's that incorporate catch by period.*/ /*Order of standard parameters for p, surv, tag, migrate, report (18 standard parameters for saturated model, report is also by tag type 1,2,3,4). Period State 2 1 2 2 3 1 3 2 4 1 4 2 5 1 5 2 6 1 6 2 7 1 7 2 8 1 8 2 9 1 9 2 10 1 10 2 */ options noovp nodate; title 'Yellowtail Flounder - migration likelihood ' ; proc iml symsize=500 worksize=300; show space; reset nocenter; * -------------- Global variables ---------------; %let model = 26; /* specifies the design matrices */ %let nstrata= 2; /* number of strata*/ %let nsample= 10; /* number of sample times*/ %let ntype=4; /* number of tag types*/ /* model Description*/ /* 1 p,phi,tag,mig,report*/ /* 2 p_s,phi_s,tag_s,mig_s,report_s*/ /* 3 p_st,phi_st,tag_st,mig_st,report_st*/ /* 4 p_s,phi_s,tag,mig_st,report_k constant tag, report varies over tag type*/ /* 5 p_s,phi_s,tag,mig_s,report_k*/ /* 6 p_st,phi_s,tag,mig_s,report*/ /* 7 p_s,phi,tag,mig_s,report_k*/ /* 8 p_s,phi_st,tag,mig_s,report_s*/ /* 9 p_s,phi,tag,mig_s,report*/ /* 10 p_s,phi_st,tag,mig_s,report_k12*/ /* 11 p_t,phi_st,tag,mig_s,report*/ /* 12 p,phi_st,tag,mig_s,report*/ /* 13 p_s,phi_st,tag,mig_s,report*/ /* 14 p_st,phi_st,tag,mig_s,report*/ /* Final model June 06*/ /* 15 p_t,phi,tag,mig_s,report exploitation varies over time, survival constant*/ /* 16 p_st, phi, tag, mig_s, report */ /* 17 p_st, phi_st, tag, mig, report*/ /* 18 p_st, phi, tag, mig, report*/ /* 19 p_st, phi, tag, mig_s, report*/ /*Same as model 16*/ /* 20 p_st, phi_t, tag, mig_s, report*/ /* 21 p_st, phi_ts2, tag, mig_s, report*/ /*survival differs by stratum in 2001 and 2004, all other years will be held constant over stratum*/ /* 22 p_st, phi_ts3, tag, mig_s, report*/ /*survival differs by stratum in 2004, all other years will be held constant over stratum*/ /* 23 p_st, phi_s, tag, mig_st, report */ /*Aug 2007 added to deterimine why model 14 has a high M estimate in 2004*/ /* 24 p_st, phi_st, tag, mig_st, report*/ /* 25 p_stp, phi_t, tag, mig_s, report*/ /*Oct 2007, exploitation varies by period*/ /* 26 p_stp, phi, tag, mig_s, report*/ /*Nov 2007, exploitation varies by period, survival constant over years*/ /* 27 p_stp, phi, tag, mig, report*/ /*Dec 2007, exploitation varies by period, survival constant*/ /* 28 p_st, phi_st, tag, mig_st, report_st*/ /*Dec 2007, */ /* 29 p_st, phi, tag, mig_st, report*/ /*Dec 2007, */ /* 30 p_st, phi, tag, mig_t, report*/ /*Dec 2007, */ /* 31 p_t, phi, tag, mig, report*/ /*Dec 2007, */ /* Additional by period models*/ /* 32 p_stp,phi_st,tag_st,mig_st,report_st*/ /* 33 p_stp,phi_st,tag,mig_st,report_st*/ /* 34 p_stp,phi_st,tag,mig_st,report*/ /* 35 p_stp,phi,tag,mig_st,report*/ /* 36 p_stp,phi,tag,mig,report*/ /* 37 p_p,phi,tag,mig,report*/ *-------------- mylogit ------------------------; * Looks at a column vector and if the values are between 0 and 1 then the logit function is evaluated; * If x=0 or x=1 (the boundary) values of -10 or 10 are forced; start mylogit(x,ans); n=nrow(x); *length of x matrix; do i=1 to n; if x[i]=0 then ans[i]=-10; else if x[i]=1 then ans[i]=10; else ans[i]=log(x[i]/(1-x[i])); *logit transform; end; finish mylogit; *-------------- qaic function ------------------------; start qaicf(R,R_expected,overd,sum_nbeta,alpha); *Calculates the QAIC value for a given model.; CALL calc_loglik(sum_loglik, R_expected,R,alpha); *overd= 1.4032493; /*overdispersion of model 32*/ overd= 1.401074; /*overdispersion from final model 26*/ *q= -2*sum_loglik/overd +2*(sum_nbeta+1); q= -2*sum_loglik +2*(sum_nbeta+1); *for NB model use AIC not QAIC; return(q); finish qaicf; *-------------- create residual function ------------------------; start Residual_calc(residual,cutoff,max_cutoff,overd,deviance,df,D,residual_strata,residual_type, residual_sample)global(sum_nbeta,N_double,std_p, std_surv, std_tag, std_report, std_migrate,R); count=0; do i=1 to &nsample-1 by 2; *released every 2 periods; do j=i+1 to &nsample; do k=1 to &nstrata; *For each recovery I have an nstrata X nstrata matrix; count=count+1; end; end; end; * initialize the expected recovery matrix to zero; R_expected=j(&ntype*count,&nstrata,0); * Initialize the matrix for the expected recoveries to zero; * calculate the expected recoveries; CALL expected_recovery(R_expected, N_double,std_p, std_surv, std_tag, std_report, std_migrate); *Calculated elementwise matrix residuals; n=nrow(R_expected); do i=1 to n; do j=1 to &nstrata; if R_expected[i,j]>0 then residual[i,j]=(R[i,j] -R_expected[i,j])/sqrt(R_expected[i,j]); else residual[i,j]=(R[i,j]-R_expected[i,j]); end; end; *Calculated elementwise matrix overdispersion; overd=0; div=0; do i=1 to n; do j=1 to &nstrata; if R_expected[i,j]>cutoff then do; if R_expected[i,j]3) then do; counter=0; do k=1 to &ntype; do m=1 to &nsample -1 by 2; do n=m+1 to &nsample; do l=1 to &nstrata; counter=counter+1; if counter=i then print 'type ' k 'release period ' m 'recover period ' n 'release stratum ' l 'recover stratum ' j 'resid ' res 'Rexpected' Rexp 'Robs' Robs; end; end; end; end; end; end; end; *print R R_expected residual overdispersion; finish Residual_calc; *-------------- transform beta to std function ------------------------; start transform_beta_to_std(beta_val,std_val); n=max(ncol(beta_val),nrow(beta_val)); *should be either a column or row vector; *std_val[1:n]=exp(beta_val[1:n])/(1+exp(beta_val[1:n])); do i=1 to n; std_val[i]=exp(max(-20,min(20,beta_val[i])))/(1+exp(max(-20,min(20,beta_val[i])))); end; finish transform_beta_to_std; *-------------- transform beta to alpha function ------------------------; start transform_beta_to_alpha(beta_val,alpha_val); * Transforms the beta value to the standard alpha overdispersion parameter. Note that the minimum beta value can be -20 as I am using a log transform so that alpha is not a negative number; alpha_val[1]=exp(max(-20,min(20,beta_val[1]))); finish transform_beta_to_alpha; *-------------- calculate loglikelihood function ------------------------; start calc_loglik(sum_loglik, R_expected, R, alpha); * Calculates the negative binomial loglikelihood, R_expected is the expected count, R is the observed recovieries, alpha is the overdisperison parameter. If expected counts are <=0, log(expected) is set to 0.; *print 'inside calc_loglik'; *print 'alpha in calc_loglik' alpha; index=1; sum_loglik=0; do tag=1 to &ntype; do rel=1 to &nsample-1 by 2; do rec=rel+1 to &nsample; do rel_strata=1 to &nstrata; do rec_strata=1 to &nstrata; *Negative binomial loglikelihood; if R_expected[index,rec_strata]>0 then loglik=R[index,rec_strata]*log(alpha*R_expected[index,rec_strata]) - (R[index,rec_strata] + 1/alpha)* log(1+alpha*R_expected[index,rec_strata]) + lgamma(R[index,rec_strata] + 1/alpha) - lgamma(R[index,rec_strata]+1) - lgamma(1/alpha); else loglik=lgamma(R[index,rec_strata] + 1/alpha) - lgamma(R[index,rec_strata]+1) - lgamma(1/alpha); *Poisson log likelihood * loglik= -R_expected[index,rec_strata]; * if R_expected[index,rec_strata]>0 then loglik=loglik+R[index,rec_strata]* * log(R_expected[index,rec_strata]); * else loglik=loglik+R[index,rec_strata]*log(1); sum_loglik=sum_loglik+loglik; end; index=index+1; end; end; end; end; finish calc_loglik; *------------------------------individual deviances-------------------------; start calc_deviance(dev, R_expected, R); * Calculates the individual poisson cell deviances, R_expected is the expected count, R is the observed recovieries. If expected counts are <=0, they are set to 1 so that the log(1)=0.; * Find the number of rows of the release matrix associated with each tag type; count=nrow(R); dev=j(count,&nstrata,0); * Initialize the matrix for the deviances to zero; index=1; do tag=1 to &ntype; do rel=1 to &nsample-1 by 2; do rec=rel+1 to &nsample; do rel_strata=1 to &nstrata; do rec_strata=1 to &nstrata; dev[index,rec_strata]= -R_expected[index,rec_strata]; if R_expected[index,rec_strata]>0 then dev[index,rec_strata]=dev[index,rec_strata]+R[index,rec_strata]* log(R_expected[index,rec_strata]); else dev[index,rec_strata]=dev[index,rec_strata]+R[index,rec_strata]*log(1); end; index=index+1; end; end; end; end; finish calc_deviance; *-------------- loglikelihood function ------------------------; start loglikelihood(pack)global(nbeta_p, nbeta_surv, nbeta_tag, nbeta_migrate, nbeta_report, nbeta_alpha,nstd_capture, nstd_surv, nstd_tag, nstd_migrate, nstd_report,nstd_alpha, p_design, surv_design, tag_design, migrate_design, report_design, alpha_design, N_double,R); * print 'inside likelihood'; * unpack betas; call unpack_beta(nbeta_p,nbeta_surv,nbeta_tag,nbeta_migrate,nbeta_report,nbeta_alpha,pack,beta_p,beta_surv,beta_tag, beta_migrate,beta_report,beta_alpha); *print beta_alpha; * transform betas to standard; new_p=beta_p; new_surv=beta_surv; new_migrate=beta_migrate; new_tag=beta_tag; new_report=beta_report; new_alpha=beta_alpha; *print 'transforming beta'; call transform_beta_to_std(beta_p,new_p); call transform_beta_to_std(beta_surv,new_surv); call transform_beta_to_std(beta_migrate,new_migrate); call transform_beta_to_std(beta_tag,new_tag); call transform_beta_to_std(beta_report,new_report); call transform_beta_to_alpha(beta_alpha,new_alpha); *print new_alpha; * multiply design matrices by my paranmeter matrices; big_p=j(nstd_capture,nbeta_p,0); big_surv=j(nstd_surv,nbeta_surv,0); big_tag=j(nstd_tag,nbeta_tag,0); big_migrate=j(nstd_migrate,nbeta_migrate,0); big_report=j(nstd_report,nbeta_report,0); big_alpha=j(nstd_alpha,nbeta_alpha,0); big_p=p_design*new_p; big_surv=surv_design*new_surv; big_tag=tag_design*new_tag; big_migrate=migrate_design*new_migrate; big_report=report_design*new_report; big_alpha=alpha_design*new_alpha; * print big_alpha; * Create 2x2xn standard parameter matrices so that I can pull of a 2x2 matric for multiplication in my expected values.; * Note that the index for p will start at 1 not at 2.; CALL create_standard_matrix(big_p,std_p); CALL create_standard_matrix(big_surv,std_surv); CALL create_standard_matrix(big_tag,std_tag); CALL create_standard_matrix_report(big_report,std_report); CALL create_standard_matrix_migrate(big_migrate,std_migrate); * Calculate expected recovery; * Find the number of rows of the release matrix associated with each tag type; count=0; do i=1 to &nsample-1 by 2; *released every 2 periods; do j=i+1 to &nsample; do k=1 to &nstrata; *For each recovery I have an nstrata X nstrata matrix; count=count+1; end; end; end; R_expected=j(&ntype*count,&nstrata,0); * Initialize the matrix for the expected recoveries to zero; *print 'into expected recovery calculation'; /*Hard code in lambda value from Poisson model = 0.988*/ /*std_report=0.988;*/ CALL expected_recovery(R_expected, N_double,std_p, std_surv, std_tag, std_report, std_migrate); *print 'calculated expected recovery'; * Calculate the log likelihood for the poisson model; * print 'alpha into calc_loglik' new_alpha; CALL calc_loglik(sum_loglik, R_expected, R,new_alpha); * print sum_loglik; * print 'finish loglikelihood'; return(sum_loglik); finish loglikelihood; *-------------- expected recovery function ------------------------; start expected_recovery(R_expected, N_double, std_p, std_surv, std_tag, std_report, std_migrate); I=diag({1 1}); *identity matrix; temp=j(&nstrata,&nstrata,0); report_other=j(&nstrata,&nstrata,0); *if assume high reward reporting rate is something other than 1, multiply by this matrix; report_other[1,1]=1; report_other[2,2]=1; * R_expected = number of releases, as must initially multiply by N; index=1; *index=1,...,360; index2=1; *index2=1,...,18; do k=1 to &ntype; do release=1 to &nsample-1 by 2; do recover=release+1 to &nsample; R_expected[index:index+&nstrata-1,1:&nstrata]=N_double[index2:index2+&nstrata-1,1:&nstrata]; index=index+&nstrata; end; index2=index2+&nstrata; end; end; index=1; *Index for R_expected and report matrices, index=1,...,200; do k=1 to &ntype; do release=1 to &nsample-1 by 2; *index3 is the index for the standard report matrix which differs by recovery period and tag type. So for tag type 2 I have to increase by nstrata*nsample-1 indeces. Not sure if i have tagtype 2 in the reporting matrix or not. tagtype release index3 1 1 1 1 3 5 1 9 17 2 1 19 2 9 35 3 1 37 3 9 53 4 1 55 4 9 71; index3=(k-1)*((&nsample-1)*&nstrata)+release*2-1; *index3=1,...,72; do recover=release+1 to &nsample; *index 2 is the index for the standard parameter matrices except for report. It is the release time multiplied by 2 (as there are 2 periods) -1 to get the correct index. So we have release index2 1 1 3 5 5 9 7 13 9 17; index2=release*2-1; *1. single or high reward tag; if k=1 | k=4 then do; *single or high reward tag; * using a temporary nstrata X nstrata matrix to do the matrix multiplications; temp[1:&nstrata, 1:&nstrata]=R_expected[index:index+&nstrata-1, 1:&nstrata]; do l=release to recover-1; * sample points between release and recover; if l=recover-1 then do; * if at recovery sample point, then multiply by capture and reporting; *single tag; if k=1 then temp=temp*std_surv[index2:index2+&nstrata-1,1:&nstrata] * std_tag[index2:index2+&nstrata-1,1:&nstrata] * std_migrate[index2:index2+&nstrata-1,1:&nstrata]* std_p[index2:index2+&nstrata-1,1:&nstrata]* std_report[index3:index3+&nstrata-1,1:&nstrata]; else temp=temp*std_surv[index2:index2+&nstrata-1,1:&nstrata]* std_tag[index2:index2+&nstrata-1,1:&nstrata]*std_migrate[index2:index2+&nstrata-1,1:&nstrata]* std_p[index2:index2+&nstrata-1,1:&nstrata]*report_other[1:&nstrata,1:&nstrata]; *high reward has report=1; end;*if; *if not at recovery sample point then multiply by (I-capture); else temp=temp*std_surv[index2:index2+&nstrata-1,1:&nstrata]* std_tag[index2:index2+&nstrata-1,1:&nstrata]*std_migrate[index2:index2+&nstrata-1,1:&nstrata]* (I-std_p[index2:index2+&nstrata-1,1:&nstrata]); index2=index2+&nstrata; end;*do; R_expected[index:index+&nstrata-1, 1:&nstrata]=temp; index=index+&nstrata; index3=index3+&nstrata; end;*if; *2. Double tag; if k=3 then do; temp[1:&nstrata, 1:&nstrata]=R_expected[index:index+&nstrata-1, 1:&nstrata]; do l=release to recover-1; * sample points between release and recover; * if at recovery sample point, then multiply by capture and reporting; if l=recover-1 then temp=temp*std_surv[index2:index2+&nstrata-1,1:&nstrata]* std_tag[index2:index2+&nstrata-1,1:&nstrata]*std_tag[index2:index2+&nstrata-1,1:&nstrata]* std_migrate[index2:index2+&nstrata-1,1:&nstrata]* std_p[index2:index2+&nstrata-1,1:&nstrata]* std_report[index3:index3+&nstrata-1,1:&nstrata]; *if not at recovery sample point then multiply by (I-capture); else temp=temp*std_surv[index2:index2+&nstrata-1,1:&nstrata]* std_tag[index2:index2+&nstrata-1,1:&nstrata]*std_tag[index2:index2+&nstrata-1,1:&nstrata]* std_migrate[index2:index2+&nstrata-1,1:&nstrata]*(I-std_p[index2:index2+&nstrata-1,1:&nstrata]); index2=index2+&nstrata; end;*do l; R_expected[index:index+&nstrata-1, 1:&nstrata]=temp; index=index+&nstrata; index3=index3+&nstrata; end;*if k=3; *3. Single tag released as double tag; if k=2 then do; temp=I; *Here temp is initialized to 1 as I will multiply by the number released at the very end; Lam_temp=j(&nstrata,&nstrata,0); *1. Possible times of tag lost; do time_tag_lost=release to recover-1; temp=I; *Here temp is initialized to 1 as I will multiply by the number released at the very end; index2=release*2-1; *2. Before time of tag loss; do j=release to time_tag_lost-1; temp=temp*std_surv[index2:index2+&nstrata-1,1:&nstrata]* std_tag[index2:index2+&nstrata-1,1:&nstrata]*std_tag[index2:index2+&nstrata-1,1:&nstrata]* std_migrate[index2:index2+&nstrata-1,1:&nstrata]*(I-std_p[index2:index2+&nstrata-1,1:&nstrata]); index2=index2+&nstrata; end;*j; *3. At time of tag loss; *Tag lost at last time period; if time_tag_lost=recover-1 then temp=temp*std_surv[index2:index2+&nstrata-1,1:&nstrata]* 2*std_tag[index2:index2+&nstrata-1,1:&nstrata]*(I-std_tag[index2:index2+&nstrata-1,1:&nstrata])* std_migrate[index2:index2+&nstrata-1,1:&nstrata]* std_p[index2:index2+&nstrata-1,1:&nstrata]* std_report[index3:index3+&nstrata-1,1:&nstrata]; *Tag lost before recovery period; else do; temp=temp*std_surv[index2:index2+&nstrata-1,1:&nstrata]* 2*std_tag[index2:index2+&nstrata-1,1:&nstrata]*(I-std_tag[index2:index2+&nstrata-1,1:&nstrata])* std_migrate[index2:index2+&nstrata-1,1:&nstrata]*(I-std_p[index2:index2+&nstrata-1,1:&nstrata]); index2=index2+&nstrata; *4. After time of tag lost but not at recovery time; do j=time_tag_lost+1 to recover-1; *5. At time of recovery; if j=(recover-1) then temp=temp*std_surv[index2:index2+&nstrata-1,1:&nstrata]* std_tag[index2:index2+&nstrata-1,1:&nstrata]*std_migrate[index2:index2+&nstrata-1,1:&nstrata]* std_p[index2:index2+&nstrata-1,1:&nstrata]*std_report[index3:index3+&nstrata-1,1:&nstrata]; if j=recover-1 then temp=temp; else temp=temp*std_surv[index2:index2+&nstrata-1,1:&nstrata]* std_tag[index2:index2+&nstrata-1,1:&nstrata]*std_migrate[index2:index2+&nstrata-1,1:&nstrata]* (I-std_p[index2:index2+&nstrata-1,1:&nstrata]); index2=index2+&nstrata; end; *do; end;*else; *6. Add up all ways of losing a tag; if i=release then Lam_temp=temp; else Lam_temp=Lam_temp+temp; end; *time_tag_lost; R_expected[index:index+&nstrata-1,1:&nstrata]=R_expected[index:index+&nstrata-1, 1:&nstrata]*Lam_temp; index=index+&nstrata; index3=index3+&nstrata; end;*if k=2; end;*recover; end;*release; end; *ntype; finish expected_recovery; *-------------- create parameter matrix function ------------------------; start create_standard_matrix(big,x); *Creates the 2*nsample X 2 matrix of standard parameter value. I will pick off individual 2x2 matrices to multiply out for the expected values.; index=1; x=j(&nsample*&nstrata-2,&nstrata,0); do until(index=&nsample*&nstrata-1); do i=1 to &nstrata; x[index,i]=big[index]; index=index+1; end; end; finish create_standard_matrix; *-------------- create migration matrix function ------------------------; start create_standard_matrix_migrate(big,x); *Creates the 2*nsample X 2 matrix of standard parameter value. I will pick off individual 2x2 matrices to multiply out for the expected values.; *Setting the off-diagonals; index=1; index2=1; x=j(&nsample*&nstrata-2,&nstrata,0); do until(index2=&nsample*&nstrata-1); do i=1 to &nstrata; do j=1 to &nstrata; if i^=j then do; x[index2,j]=big[index]; index=index+1; end; end; index2=index2+1; end; end; * Setting the diagonals; index=1; do until(index=&nsample*&nstrata-1); do i= 1 to &nstrata; x[index,i]=1-sum(x[index,1:&nstrata]); index=index+1; end; end; finish create_standard_matrix_migrate; *-------------- create report matrix function ------------------------; start create_standard_matrix_report(big,x); *Creates the ntype*(2*nsample) X 2 matrix of standard parameter value. I will pick off individual 2x2 matrices to multiply out for the expected values.; index=1; x=j(&ntype*(&nsample*&nstrata-2),&nstrata,0); do until(index=&ntype*(&nsample*&nstrata-2)+1); do i=1 to &nstrata; x[index,i]=big[index]; index=index+1; end; end; finish create_standard_matrix_report; *-------------- standard to beta function ------------------------; *Find inital values for the betas using least squares; *Given x=standard parameters, design=design matrix, using a logistic link ie y=logistic(phi)=XB; *Then beta=(x`x)^-1 x`y which is from standard least squares regression.; start std_to_beta(x,design, beta); y=x;/*initialize y (transformed standard parameters) to something, same length as x*/ CALL mylogit(x,y); beta=inv(design`*design)*design`*y; *could also use beta=ginv(design)*y; finish std_to_beta; *-------------- standard to beta function ------------------------; *Find inital values for the betas using least squares; *Given x=standard parameters, design=design matrix, using a log link ie y=log(phi)=XB; *Then beta=(x`x)^-1 x`y which is from standard least squares regression.; *Here I only have one alpha which is the overdispersion parameter so this may be unnecessary.; start alpha_to_beta(x,design, beta); y=x;/*initialize y (transformed standard parameters) to something, same length as x*/ y=log(x); beta=inv(design`*design)*design`*y; *could also use beta=ginv(design)*y; finish alpha_to_beta; *-------------- transform beta vector to std vector ------------------------; start transform_to_std(X,beta); * Transform the beta estimates to standard estimates by first multiplying by the large design matrix; all_beta=X*beta`; std_parms=all_beta; CALL transform_beta_to_std(all_beta,std_parms); return(std_parms); finish transform_to_std; *-------------- beta pack -----------------------; start beta_pack(nbeta_p,nbeta_surv,nbeta_tag,nbeta_migrate,nbeta_report,nbeta_alpha,beta_p,beta_surv,beta_tag,beta_migrate,beta_report,beta_alpha,pack); pack[1:nbeta_p]=beta_p; pack[nbeta_p+1:nbeta_p+nbeta_surv]=beta_surv; pack[nbeta_p+nbeta_surv+1:nbeta_p+nbeta_surv+nbeta_tag]=beta_tag; pack[nbeta_p+nbeta_surv+nbeta_tag+1:nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate]=beta_migrate; pack[nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+1:nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+nbeta_report]=beta_report; pack[nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+nbeta_report+1:nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+nbeta_report+nbeta_alpha]=beta_alpha; finish beta_pack; *-------------- beta unpack -----------------------; start unpack_beta(nbeta_p,nbeta_surv,nbeta_tag,nbeta_migrate,nbeta_report,nbeta_alpha,pack,beta_p,beta_surv,beta_tag,beta_migrate,beta_report,beta_alpha); beta_p=pack[1:nbeta_p]; beta_surv=pack[nbeta_p+1:nbeta_p+nbeta_surv]; beta_tag =pack[nbeta_p+nbeta_surv+1:nbeta_p+nbeta_surv+nbeta_tag]; beta_migrate=pack[nbeta_p+nbeta_surv+nbeta_tag+1:nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate]; beta_report=pack[nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+1:nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+nbeta_report]; beta_alpha=pack[nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+nbeta_report+1:nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+nbeta_report+nbeta_alpha]; finish unpack_beta; *-------------- create big design matrix -----------------------; start make_large_design(p_design, surv_design, tag_design, migrate_design, report_design,alpha_design, big_design) global(nbeta_p, nbeta_surv, nbeta_tag, nbeta_migrate, nbeta_report, nstd_capture, nstd_surv, nstd_tag, nstd_migrate, nstd_report, sum_nbeta, sum_nstd); big_design[1:nstd_capture,1:nbeta_p]=p_design; big_design[(nstd_capture+1):(nstd_capture+nstd_surv), (nbeta_p+1):(nbeta_p+nbeta_surv)]=surv_design; big_design[(nstd_capture+nstd_surv+1):(nstd_capture+nstd_surv+nstd_tag), (nbeta_p+nbeta_surv+1):(nbeta_p+nbeta_surv+nbeta_tag)]=tag_design; big_design[(nstd_capture+nstd_surv+nstd_tag+1):(nstd_capture+nstd_surv+nstd_tag+nstd_migrate), (nbeta_p+nbeta_surv+nbeta_tag+1):(nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate)]=migrate_design; big_design[(nstd_capture+nstd_surv+nstd_tag+nstd_migrate+1):(nstd_capture+nstd_surv+nstd_tag+nstd_migrate+ nstd_report), (nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+1):(nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+nbeta_report)]=report_design; big_design[(nstd_capture+nstd_surv+nstd_tag+nstd_migrate+nstd_report+1):sum_nstd, (nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+nbeta_report+1):sum_nbeta]=alpha_design; finish make_large_design; *-------------- create transformed parameters -----------------------; start transformed_parameters(F, M, tau, stderr_TP)global(new_p, new_surv, new_tag, stderr, nstd_capture, nstd_surv, nstd_tag); F=-log(1-2*new_p); M=-2*log(new_surv); tau=-2*log(new_tag); stderr_TP[1:nstd_capture]= (2/(1-2*new_p[1:nstd_capture]))*stderr[1:nstd_capture]; stderr_TP[(nstd_capture+1):(nstd_capture+nstd_surv+nstd_tag)]= (2/new_surv[(nstd_capture+nstd_surv+nstd_tag)])*stderr[(nstd_capture+1):(nstd_capture+nstd_surv+nstd_tag)]; finish transformed_parameters; *-------------- begin program ------------------------------------------------------------------------------------------; if &model=1 then do; print 'Using model 1: p, phi, tag, Q, report'; * Create model 1 (all parameters constant)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {capture}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_c"; do data; input capture; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {migrate}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_c"; do data; input migrate; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {surv}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input surv; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=2 then do; print 'Using model 2: p_s, phi_s, tag_s, Q_s, report_s'; * Create model 2 (all parameters vary by stratum)- read in design matrices; create design_tag var {t1 t2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_s"; do data; input t1 t2; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {c1 c2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_s"; do data; input c1 c2; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{r1 r2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_s"; do data; input r1 r2; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_s"; do data; input s1 s2; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=3 then do; print 'Using model 3: p_st, phi_st, tag_st, Q_st, report_st'; * Create model 3 (all parameters vary by stratum and time)- read in design matrices; create design_tag var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2 m3 m4 m5 m6 m7 m8 m9 m10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_st"; do data; input m1 m2 m3 m4 m5 m6 m7 m8 m9 m10; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model =4 then do; print 'Using model 4: p_s, phi_s, tag, Q_st, report_k'; * Create model 4 (p_s, phi_s, tag, Q_st, report_k)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {c1 c2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_s"; do data; input c1 c2; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{r1 r2 r3 r4}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_k"; do data; input r1 r2 r3 r4; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2 m3 m4 m5 m6 m7 m8 m9 m10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_st"; do data; input m1 m2 m3 m4 m5 m6 m7 m8 m9 m10; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_s"; do data; input s1 s2; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=5 then do; print 'Using model 5: p_s, phi_s, tag, Q_s, report_k'; * Create model 5 (p_s, phi_s, tag, Q_s, report_k)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {c1 c2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_s"; do data; input c1 c2; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{r1 r2 r3 r4}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_k"; do data; input r1 r2 r3 r4; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_s"; do data; input s1 s2; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=6 then do; print 'Using model 6: p_st, phi_s, tag, Q_s, report'; * Create model 6 (p_st, phi_s, tag, Q_s, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {c1 c2 c3 c4 c5 c6 c7 c8 c9 c10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input c1 c2 c3 c4 c5 c6 c7 c8 c9 c10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_s"; do data; input s1 s2; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=7 then do; print 'Using model 7: p_s, phi, tag, Q_s, report_k'; * Create model 1 (p_s, phi, tag, Q_s, report_k)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {c1 c2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_s"; do data; input c1 c2; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{r1 r2 r3 r4}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_k"; do data; input r1 r2 r3 r4; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {surv}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input surv; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=8 then do; print 'Using model 8: p_s, phi_st, tag, Q_s, report_s'; * Create model 8 (p_s, phi_st, tag, Q_s, report_s)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {c1 c2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_s"; do data; input c1 c2; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{r1 r2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_s"; do data; input r1 r2; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=9 then do; print 'Using model 9: p_s, phi, tag, Q_s, report'; * Create model 9 (p_s, phi, tag, Q_s, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {c1 c2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_s"; do data; input c1 c2; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {surv}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input surv; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=10 then do; print 'Using model 10: p_s, phi_st, tag, Q_s, report_k12'; * Create model 1 (p_s, phi_st, tag, Q_s, report_k12)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {c1 c2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_s"; do data; input c1 c2; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{r1 r2 r3}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_k12"; do data; input r1 r2 r3; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=11 then do; print 'Using model 11: p_t, phi_st, tag, Q_s, report'; * Create model 11 (p_t, phi_st, tag, Q_s, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {c1 c2 c3 c4 c5}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_t"; do data; input c1 c2 c3 c4 c5; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=12 then do; print 'Using model 12: p, phi_st, tag, Q_s, report'; * Create model 12 (p, phi_st, tag, Q_s, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {capture}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_c"; do data; input capture; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=13 then do; print 'Using model 13: p_s, phi_st, tag, Q_s, report'; * Create model 13 (p_s, phi_st, tag, Q_s, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {c1 c2 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_s"; do data; input c1 c2; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=14 then do; print 'Using model 14: p_st, phi_st, tag, Q_s, report'; * Create model 14 (p_st, phi_st, tag, Q_s, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=15 then do; print 'Using model 15: p_t, phi, tag, Q_s, report'; * Create model 15 (p_t, phi, tag, Q_s, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {c1 c2 c3 c4 c5}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_t"; do data; input c1 c2 c3 c4 c5; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {surv}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input surv; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=16 then do; print 'Using model 16: p_t, phi, tag, Q_s, report'; * Create model 16 (p_t, phi, tag, Q_s, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {surv}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input surv; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=17 then do; print 'Using model 17: p_st, phi_st, tag, Q, report'; * Create model 17 (p_st, phi_st, tag, Q, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_c"; do data; input m; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=18 then do; print 'Using model 18: p_st, phi, tag, Q, report'; * Create model 18 (p_st, phi, tag, Q, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_c"; do data; input m; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input s; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=19 then do; print 'Using model 18: p_st, phi, tag, Q, report'; * Create model 18 (p_st, phi, tag, Q, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input s; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=20 then do; print 'Using model 20: p_st, phi, tag, Q, report'; * Create model 20 (p_st, phi, tag, Q, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_t"; do data; input s1 s2 s3 s4 s5; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=21 then do; print 'Using model 20: p_st, phi, tag, Q, report'; * Create model 21 (p_st, phi, tag, Q, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6 s7}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_ts2"; do data; input s1 s2 s3 s4 s5 s6 s7; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=22 then do; print 'Using model 22: p_st, phi_ts3, tag, Q, report'; * Create model 22 (p_st, phi_ts3, tag, Q, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_ts3"; do data; input s1 s2 s3 s4 s5 s6; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; /* 23 p_st, phi_s, tag, mig_st, report */ /*Aug 2007 added to deterimine why model 14 has a high M estimate in 2004*/ else if &model=23 then do; print 'Using model 23: p_st, phi_ts3, tag, Q, report'; * Create model 23 (p_st, phi_s, tag, Q_st, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2 m3 m4 m5 m6 m7 m8 m9 m10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_st"; do data; input m1 m2 m3 m4 m5 m6 m7 m8 m9 m10; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_s"; do data; input s1 s2 ; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; /* 24 p_st, phi_st, tag, mig_st, report*/ else if &model=24 then do; print 'Using model 24: p_st, phi_st, tag, Q_st, report'; * Create model 24 (p_st, phi_st, tag, Q_st, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2 m3 m4 m5 m6 m7 m8 m9 m10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_st"; do data; input m1 m2 m3 m4 m5 m6 m7 m8 m9 m10; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_surv; use design_surv; READ all INTO surv_design; END; else if &model=25 then do; print 'Using model 25: p_st, phi_t, tag, Q_s, report'; * Create model 25 (p_st, phi_st, tag, Q_s, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_stp"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_t"; do data; input s1 s2 s3 s4 s5; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=26 then do; print 'Using model 26: p_st, phi_t, tag, Q_s, report'; * Create model 26 (p_st, phi, tag, Q_s, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_stp"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_s"; do data; input m1 m2; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input s; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=27 then do; print 'Using model 27: p_stp, phi, tag, Q, report'; * Create model 27 (p_stp, phi, tag, Q, report)- read in design matrices; create design_tag var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input tag; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_stp"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input report; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_c"; do data; input m1; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input s; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=28 then do; print 'Using model 28: p_st, phi_st, tag, Q_st, report_st'; * Create model 28 (p_stp, phi, tag, Q, report)- read in design matrices; create design_tag var {t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input t1; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2 m3 m4 m5 m6 m7 m8 m9 m10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_st"; do data; input m1 m2 m3 m4 m5 m6 m7 m8 m9 m10; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; READ all INTO migrate_design; create design_surv var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=29 then do; print 'Using model 29: p_st, phi, tag, Q_st, report'; * Create model 29 (p_st, phi, tag, Q, report)- read in design matrices; create design_tag var {t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input t1; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input t1; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2 m3 m4 m5 m6 m7 m8 m9 m10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_st"; do data; input m1 m2 m3 m4 m5 m6 m7 m8 m9 m10; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; READ all INTO migrate_design; create design_surv var {s1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input s1; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=30 then do; print 'Using model 30: p_st, phi, tag, Q_t, report'; * Create model 30 (p_st, phi, tag, Q_t, report)- read in design matrices; create design_tag var {t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input t1; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input t1; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {m1 m2 m3 m4 m5}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_t"; do data; input m1 m2 m3 m4 m5; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; READ all INTO migrate_design; create design_surv var {s1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input s1; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=31 then do; print 'Using model 31: p_st, phi, tag, Q, report'; * Create model 31 (p_st, phi, tag, Q, report)- read in design matrices; create design_tag var {t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input t1; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {t1 t2 t3 t4 t5}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_t"; do data; input t1 t2 t3 t4 t5; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input t1; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {migrate}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_c"; do data; input migrate; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {s1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input s1; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=32 then do; print 'Using model 32: p_sp, phi_st, tag_st, Q__st, report_st'; * Create model 32 (p_sp, phi_st, tag_st, Q_st, report_st)- read in design matrices; create design_tag var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_stp"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=33 then do; print 'Using model 33: p_sp, phi_st, tag_st, Q__st, report_st'; * Create model 33 (p_sp, phi_st, tag_st, Q_st, report_st)- read in design matrices; create design_tag var {t1 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input t1; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_stp"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=33 then do; print 'Using model 33: p_sp, phi_st, tag_st, Q__st, report_st'; * Create model 33 (p_sp, phi_st, tag_st, Q_st, report_st)- read in design matrices; create design_tag var {t1 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input t1; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_stp"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=34 then do; print 'Using model 34: p_sp, phi_st, tag_st, Q__st, report'; * Create model 34 (p_sp, phi_st, tag_st, Q_st, report)- read in design matrices; create design_tag var {t1 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input t1; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_stp"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input t1; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=35 then do; print 'Using model 35: p_sp, phi, tag_st, Q__st, report'; * Create model 35 (p_sp, phi, tag_st, Q_st, report)- read in design matrices; create design_tag var {t1 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input t1; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_stp"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input t1; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {t1 t2 t3 t4 t5 t6 t7 t8 t9 t10}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_st"; do data; input t1 t2 t3 t4 t5 t6 t7 t8 t9 t10; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input t1; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=36 then do; print 'Using model 36: p_sp, phi, tag, Q, report'; * Create model 36 (p_sp, phi, tag, Q, report)- read in design matrices; create design_tag var {t1 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input t1; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_stp"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input t1; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_c"; do data; input t1; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input t1; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; else if &model=37 then do; print 'Using model 37: p_sp, phi, tag, Q_t, report'; * Create model 37 (p_sp, phi, tag, Q_t, report)- read in design matrices; create design_tag var {t1 }; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_tag_c"; do data; input t1; append; end; close design_tag; use design_tag; READ all INTO tag_design; create design_capture var {s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_p_stp"; do data; input s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16 s17 s18; append; end; close design_capture; use design_capture; READ all INTO p_design; create design_report var{t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_report_c"; do data; input t1; append; end; close design_report; use design_report; READ all INTO report_design; create design_migrate var {t1 t2 t3 t4 t5}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_migrate_t"; do data; input t1 t2 t3 t4 t5; append; end; close design_migrate; use design_migrate; READ all INTO migrate_design; create design_surv var {t1}; infile "D:\contract\Walsh\contract2\Dec08\R Code\design_surv_c"; do data; input t1; append; end; close design_surv; use design_surv; READ all INTO surv_design; end; print 'starting alpha_design'; alpha_design=j(1,1,1); print alpha_design; * Alpha overdispersion parameter for NB ; *--- Number of standard parameters ----; nstd_tag=nrow(tag_design); nstd_capture=nrow(p_design); nstd_report=nrow(report_design); nstd_migrate=nrow(migrate_design); nstd_surv=nrow(surv_design); nstd_alpha=nrow(alpha_design); print nstd_alpha; *--- Number of beta parameters ----; nbeta_tag=ncol(tag_design); nbeta_p=ncol(p_design); nbeta_report=ncol(report_design); nbeta_migrate=ncol(migrate_design); nbeta_surv=ncol(surv_design); nbeta_alpha=1; print nbeta_tag nbeta_p nbeta_report nbeta_migrate nbeta_surv nbeta_alpha; *--- number of fish released ---; * Input file is N_11 0 ---; * 0 N_22 ---; * For Tag 1 period 1,3,5,7,9 Tag 3 period 1,3,5,7,9 Tag 4 period 1,3,5,7,9 ---; create num_rel var {N1 N2}; infile "D:\contract\Walsh\contract2\Dec08\SAS Code\num_rel.txt"; do data; input N1 N2; append; end; close num_rel; use num_rel; * Do not need num_rel=j((&ntype-1)*(&nsample/2)*&nstrata,2,100); /*num_rel={10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000, 10000 0,0 10000}; N=num_rel; */ READ all var{N1 N2}INTO N ; * This N (number of releases) does not have type 2 tags (tagged as double, recaptured as single); * This causes problems when calculating the expectations. I will lengthen N to include a spot for type 2 tags, by repeating type 3 (double tags); * So N has the following matrices: [type 1] and N_double has [type 1] [type 3] [type 2] [type 4] [type 3] [type 4]; count=0; *Number of rows for release matrix associated with each type of tag; do i=1 to &nsample-1 by 2; *Releases occur every second period, ie 1,3,5,7,9; do k=1 to &nstrata; *For each release I have an nstrata X nstrata matrix; count=count+1; end; end; nrelease_by_type=count; N_double=j(&ntype*count,&nstrata,0); *This is the matrix that includes type 2 releases (equal to type 3); index=1; do i=1 to count; *type 1, single tagged; N_double[i,1:&nstrata]=N[i,1:&nstrata]; index=index+1; end; index2=index; do i=1 to count;*type 2 double tagged, single recap; N_double[index2,1:&nstrata]=N[index2,1:&nstrata]; index2=index2+1; end; do i=1 to count; *type 3 double tagged; N_double[index2,1:&nstrata]=N[index,1:&nstrata]; index2=index2+1; index=index+1; end; do i=1 to count; *type 4 high reward; N_double[index2,1:&nstrata]=N[index,1:&nstrata]; index2=index2+1; index=index+1; end; *--- number of fish recovered ---; * Input file is R_11 R_12 ---released in stratum 1, recovered in stratum1 released in stratum 1 recovered in stratum 2; * R_21 R_22 ---released in stratum 2, recovered in stratum1 released in stratum 2 recovered in stratum 2; * For Tag 1 released period 1, recovered period 2,3,4,5,6,7,8,9,10 released period 3, recovered period 4,5,6,7,8,9,10 ... Tag 2 released period 1, recovered period 2,3,4,5,6,7,8,9,10 released period 3, recovered period 4,5,6,7,8,9,10 ... Tag 3 released period 1, recovered period 2,3,4,5,6,7,8,9,10 released period 3, recovered period 4,5,6,7,8,9,10 ... Tag 4 released period 1, recovered period 2,3,4,5,6,7,8,9,10 released period 3, recovered period 4,5,6,7,8,9,10 ... ---; create num_rec var {R1 R2}; infile "D:\contract\Walsh\contract2\Dec08\SAS Code\num_rec.txt"; do data; input R1 R2; append; end; close num_rec; use num_rec; READ all var{R1 R2}INTO R ; Rcheck=sum(R); Ncheck=sum(N_double); *print Rcheck Ncheck; *R=j(&ntype*&nstrata*&nstrata*45,&nstrata,200); * Find the number of rows of the release matrix associated with each tag type; * Only need to do this for 1 tag type; count=0; do i=1 to &nsample-1 by 2; *released every 2 periods; do j=i+1 to &nsample; do k=1 to &nstrata; *For each recovery I have an nstrata X nstrata matrix; count=count+1; end; end; end; *R_expected=j(&ntype*count,&nstrata,0); * Initialize the matrix for the expected recoveries to zero; *nrecover_by_type=count; *--- Initialize standard parameters (hard coded in) ---; create taginit var {tag}; infile "D:\contract\Walsh\contract2\Dec08\R Code\tag_initial"; do data; input tag; append; end; close taginit; use taginit; READ all INTO tag; create pinit var {p}; infile "D:\contract\Walsh\contract2\Dec08\R Code\p_initial"; do data; input p; append; end; close pinit; use pinit; READ all INTO p; create survinit var {surv}; infile "D:\contract\Walsh\contract2\Dec08\R Code\surv_initial"; do data; input surv; append; end; close survinit; use survinit; READ all INTO surv; create miginit var {migrate}; infile "D:\contract\Walsh\contract2\Dec08\R Code\migrate_initial"; do data; input migrate; append; end; close miginit; use miginit; READ all INTO migrate; create reportinit var {report}; infile "D:\contract\Walsh\contract2\Dec08\R Code\report_initial"; do data; input report; append; end; close reportinit; use reportinit; READ all INTO report; alpha=j(1,1,2); *initial alpha value; print alpha; * Initialize the beta parameters to 0; * 1 x nbeta; *print nbeta_migrate; beta_migrate=j(nbeta_migrate,1,0); beta_p=j(nbeta_p,1,0); beta_surv=j(nbeta_surv,1,0); beta_tag=j(nbeta_tag,1,0); beta_report=j(nbeta_report,1,0); beta_alpha=j(nbeta_alpha,1,0); * Initialize beta parameters by transforming the standard parameters to beta parameters; CALL std_to_beta(migrate, migrate_design, beta_migrate); CALL std_to_beta(p, p_design, beta_p); CALL std_to_beta(surv, surv_design, beta_surv); CALL std_to_beta(tag, tag_design, beta_tag); CALL std_to_beta(report, report_design, beta_report); CALL alpha_to_beta(alpha,alpha_design, beta_alpha); print beta_alpha; *Pack betas into one vector; pack=j(nbeta_p+nbeta_surv+nbeta_tag+nbeta_migrate+nbeta_report+nbeta_alpha,1,0); CALL beta_pack(nbeta_p,nbeta_surv,nbeta_tag,nbeta_migrate,nbeta_report,nbeta_alpha,beta_p,beta_surv,beta_tag,beta_migrate,beta_report,beta_alpha,pack); print 'before entering maximization'; * rc is the termination record, if greater than 0, maximization succeeded. If less than 0, something wrong.; * xres is the parameter vector that was maximized; /*This works but requires the gradient*/ optn={1 5}; /*optn[1]=1 maximization, optn[2]=2 printing options*/ *Keep the betas between -20 and 20 so that when I exponentiate the betas I don't get exp(310) which will give me overflow errors; if &model=1 then con={-20 -20 -20 -20 -20, 20 20 20 20 20}; else if &model=2 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20}; else if &model=3 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*50 parameters*/ else if &model=4 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=5 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20}; else if &model=6 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=7 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20}; else if &model=8 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=9 then con={-20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20}; else if &model=10 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=11 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=12 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=13 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 }; else if &model=14 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; /*Constrain last survival parameter =1 (beta=20) to see if confounded with last p*/ /*else if &model=14 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 20 20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; */ else if &model=15 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20}; else if &model=16 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=17 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=18 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=19 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=20 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=21 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20}; else if &model=22 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*20 parameters*/ else if &model=23 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*24 parameters*/ /*constrain last mig_2 to 0.05 in 2004*/ /*else if &model=23 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -2.94 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 -2.94 20}; */ else if &model=24 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*32 parameters*/ /*Constrain p1_2004 to 0.01*/ /*else if &model=24 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -4.59511985 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 -4.59511985 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*32 parameters*/ /*Constrain p2_2004 to 0.01*/ /*else if &model=24 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -4.59511985 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 -4.59511985 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*32 parameters*/ /*Constrain p2_2004 to 0.037*/ /*else if &model=24 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -3.259135499 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 -3.259135499 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*32 parameters*/ else if &model=25 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*27 parameters*/ else if &model=26 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*23 parameters +1 alpha*/ else if &model=27 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*22 parameters*/ else if &model=28 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*41 parameters*/ else if &model=29 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 , 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 };/*23 parameters*/ else if &model=30 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*18 parameters*/ else if &model=31 then con={ -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20};/*9 parameters*/ else if &model=32 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*58 parameters*/ else if &model=33 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*49 parameters*/ else if &model=34 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*40 parameters*/ else if &model=35 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*31 parameters*/ else if &model=36 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*22 parameters*/ else if &model=37 then con={-20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20 -20, 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20};/*26 parameters*/ print 'entering maximization'; CALL nlpcg(rc, xres, "loglikelihood",pack,optn,con); print "made it to here"; *print 'Beta parameter estimates' xres; call unpack_beta(nbeta_p,nbeta_surv,nbeta_tag,nbeta_migrate,nbeta_report, nbeta_alpha,xres,beta_p,beta_surv,beta_tag, beta_migrate,beta_report, beta_alpha); * transform betas to standard; new_p=beta_p; new_surv=beta_surv; new_migrate=beta_migrate; new_tag=beta_tag; new_report=beta_report; new_alpha=beta_alpha; call transform_beta_to_std(beta_p,new_p); call transform_beta_to_std(beta_surv,new_surv); call transform_beta_to_std(beta_migrate,new_migrate); call transform_beta_to_std(beta_tag,new_tag); call transform_beta_to_std(beta_report,new_report); call transform_beta_to_std(beta_alpha, new_alpha); print 'standard parameter estimates' new_p new_surv new_migrate new_tag new_report new_alpha; * This function gets at the hessian which can be used to get standard errors; * If don't include par, will do quick, but less precise forward differencing to get at first and second derivatives. * par is a vector of length 3, par[1] has default 1 (I left this as it is), par[2]=ij where it determines what type of * approximation is used. If j=0 then it uses forward differencing, otherwise it uses central differencing. If i=1, 2, or 3 * then the finite differences only have to do with the number of accruate digits (ie I didn't worry about this and set it to 4). * So par[2]=41. par[3] has to do with the number of accurate digits and if missing gets put to the default. I left it as missing.; par={1 41 .}; CALL nlpfdd(func,gr,hess,"loglikelihood", xres,par); /*seems to need gradient to get right variance*/ * Make a large design matrix consiting of all design matrices for all parameters; * I need this to find the covariances of the standard parameters; sum_nbeta=sum(nbeta_p,nbeta_surv,nbeta_tag,nbeta_migrate,nbeta_report,nbeta_alpha); sum_nstd=sum(nstd_capture,nstd_surv,nstd_tag,nstd_migrate,nstd_report,nstd_alpha); design=j(sum_nstd,sum_nbeta,0); CALL make_large_design(p_design, surv_design, tag_design, migrate_design, report_design,alpha_design, design); * Calculate the Covariance(XB)=X'Cov(Beta)X where Cov(beta) is the hessian from new.beta; *print hess; inv_hessian=-ginv(hess); *print inv_hessian; cov_XB=design*inv_hessian*design`; *print cov_XB; * Calculate and print correlation; *d=ginv(sqrt(diag(inv_hessian))); *corr=d*inv_hessian*d; *cname={p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 m1 m2 t r}; *print corr [COLNAME=cname rowname=cname]; *print 'correlation matrix' corr; * Make long parameter vector for the standard parameters; std_parms=transform_to_std(design, xres); * Induce the delta method to get Covariance of standard parameters; dg_dXB=std_parms#(1-std_parms); * # multiply corresponding elements of the same matrix; cov_std=diag(dg_dXB) * cov_XB * diag(dg_dXB); *print 'covariance of standard estimates'; *print cov_std; stderr=j(nrow(cov_std),1,0); n=nrow(cov_std); do i=1 to n; stderr[i]=cov_std[i,i]; if cov_std[i,i]<0 then stderr[i]=0;/*For some reason getting negative covariance estimates due to the hessian*/ else stderr[i]=sqrt(stderr[i]); end; print 'standard error estimates of the standard parameters' stderr; *------------ Create residual plot-------------------------; beta_final=xres; *final beta parameter estimates; call unpack_beta(nbeta_p,nbeta_surv,nbeta_tag,nbeta_migrate,nbeta_report, nbeta_alpha,beta_final,beta_p,beta_surv,beta_tag, beta_migrate,beta_report, beta_alpha); * Transform betas to standard parameters. This gives me p*, surv*, tag*, migrate*; * where the * indicates that there may not be as many standard parameters as for the; * saturated model. For model 1, there will be only 1 standard parameter for each component.; new_p=beta_p; new_surv=beta_surv; new_migrate=beta_migrate; new_tag=beta_tag; new_report=beta_report; new_alpha=beta_alpha; call transform_beta_to_std(beta_p,new_p); call transform_beta_to_std(beta_surv,new_surv); call transform_beta_to_std(beta_migrate,new_migrate); call transform_beta_to_std(beta_tag,new_tag); call transform_beta_to_std(beta_report,new_report); call transform_beta_to_std(beta_alpha, new_alpha); * multiply design matrices by my paraneter matrices; big_p=j(nstd_capture,nbeta_p,0); big_surv=j(nstd_surv,nbeta_surv,0); big_tag=j(nstd_tag,nbeta_tag,0); big_migrate=j(nstd_migrate,nbeta_migrate,0); big_report=j(nstd_report,nbeta_report,0); big_alpha=j(nstd_alpha,nbeta_alpha,0); big_p=p_design*new_p; big_surv=surv_design*new_surv; big_tag=tag_design*new_tag; big_migrate=migrate_design*new_migrate; big_report=report_design*new_report; big_alpha=alpha_design*new_alpha; * Create 2x2xn standard parameter matrices so that I can pull off a 2x2 matric for multiplication in my expected values.; * Note that the index for p will start at 1 not at 2.; CALL create_standard_matrix(big_p,std_p); CALL create_standard_matrix(big_surv,std_surv); CALL create_standard_matrix(big_tag,std_tag); CALL create_standard_matrix_report(big_report,std_report); CALL create_standard_matrix_migrate(big_migrate,std_migrate); CALL create_standard_matrix(big_alpha,std_alpha); *check if this will work!; * Calculate residual (o-e)/sqrt(e) and Overdisperison [(o-e)**2/e]/ df =number of cells-#estimated parameters; * Find the number of rows of the release matrix associated with each tag type; residual=j(&ntype*count,&nstrata,0); residual_strata=j(&nstrata,1,0); residual_type=j(&ntype,1,0); residual_sample=j(&nsample,1,0); CALL Residual_calc(residual,0,100,overd, deviance,df,D,residual_strata,residual_type, residual_sample); *print residual_strata residual_type residual_sample; predicted=j(&ntype*count,&nstrata,0); * Initialize the matrix for the expected recoveries to zero; * calculate the expected recoveries; CALL expected_recovery(predicted, N_double,std_p, std_surv, std_tag, std_report, std_migrate); *------------ Create individual deviances-------------------------; *Create the individual deviances and store into a file called "individ_dev". Can change the filename for each run depending on the model number. This allows me to read in the deviances from 2 different models to find out where the major differences are coming from. In particular I will do this with the constant model (1) and the changing over stratum and time model (17) as well as look at models 21 and 22 compared with 17.; dev=j(nrow(R),&nstrata,0); * Initialize the matrix for the deviances to zero; CALL calc_deviance(dev,predicted, R); /*-predicted +observed ln(predicted)*/ dev=2*dev; dev_obs=j(&ntype*count,&nstrata,0); R2=R; CALL calc_deviance(dev_obs, R2,R); /*-observed +observed ln(observed)*/ dev_obs=2*dev_obs; dev=dev_obs-dev; /*poisson deviance=2[-(observed - predicted)+ observed*ln(observed/predicted)]*/ *print dev; /* * Turn this on only when want to print out individual deviances to a file; filename out 'D:\contract\Walsh\contract2\Dec08\SAS Code\individ_dev_model14.dat'; file out; do i=1 to nrow(dev); do j=1 to ncol(dev); put(dev[i,j]) 12.9 +2 @; *'12.9'=12 columns per value with 9 decimal spots, '+2'= move over 2 columns, '@'=hold the line to write the next value; end; put; * finishes inner do loop and puts nothing so that skips to next row of the output file.; end; closefile out; */ * Outputing the indidivual deviances as well as the tag type and period.; /* filename out 'D:\contract\Walsh\contract2\Dec08\SAS Code\individ_dev.dat'; file out; index=0; do tag=1 to &ntype; *tags; do rel=1 to &nsample-1 by 2; *release period; do rec=rel+1 to &nsample; *recovery period; do rel_strata=1 to &nstrata; index=index+1; do rec_strata=1 to &nstrata; put(dev[index,rec_strata]) 12.9 +2 @; *'12.9'=12 columns per value with 9 decimal spots, '+2'= move over 2 columns, '@'=hold the line to write the next value; put(predicted[index,rec_strata]) 12.9 +2 @; put(R[index,rec_strata]) 12.9 +2 @; end; put(tag)+2 @; put(rel)+1 @; put(rec)+1 @; put(rel_strata)+1 @; put; * finishes inner do loop and puts nothing so that skips to next row of the output file.; end; end; end; end; closefile out; */ * Adjust standard errors by overdispersion; stderr_adj=sqrt(overd)*stderr; print 'Adjusted standard errors' stderr_adj 'over dispersion' overd; * Calculated the QAIC values; qaic=qaicf(R,predicted,overd,sum_nbeta); print qaic; * Need to vectorize the matrices for plotting purposes; r=j(&ntype*count*&nstrata,1,0); * Residuals; r_adj=j(&ntype*count*&nstrata,1,0); * Residuals adjusted for overdispersion; res_adj=j(&ntype*count,&nstrata,0); p=j(&ntype*count*&nstrata,1,0); * Expected count; ob=j(&ntype*count*&nstrata,1,0); * Observed count; index=1; do i=1 to &ntype*count; do j=1 to &nstrata; r[index]=residual[i,j]; r_adj[index]=residual[i,j]/sqrt(overd); res_adj[i,j]=r_adj[index]; p[index]=predicted[i,j]; ob[index]=R2[i,j]; index=index+1; end; end; *print r_adj; *print res_adj; *print predicted R2; strata=j(&nstrata,1,0); do i=1 to &nstrata; strata[i,1]=i; end; * Create the residual plots; create recov var{r_adj p }; append; create gof var{D df}; append; create strata var{residual_strata strata}; append; * Outputs the observed counts, residuals, adjusted residuals and the predicted values to a file; filename out 'D:\contract\Walsh\contract2\Dec08\SAS Code\residuals_model.dat'; file out; do i=1 to nrow(r); put(r[i]) 15.9 +2(r_adj[i]) 15.9 +2(p[i]) 15.9 +2(ob[i]) 10.0 ; *'12.9'=12 columns per value with 9 decimal spots, '+2'= move over 2 columns, '@'=hold the line to write the next value; end; put; * finishes inner do loop and puts nothing so that skips to next row of the output file.; closefile out; * Print out large residuals; n=nrow(r_adj); do i=1 to n; res=r_adj[i]; Rexp=p[i]; Robs=ob[i]; if abs(res>3) then do; counter=0; do k=1 to &ntype; do m=1 to &nsample -1 by 2; do n=m+1 to &nsample; do l=1 to &nstrata; do j=1 to &nstrata; counter=counter+1; if counter=i then print 'type ' k 'release period ' m 'recover period ' n 'release stratum ' l 'recover stratum ' j 'resid ' res 'Rexpected' Rexp 'Robs' Robs; end; end; end; end; end; end; end; proc plot data=recov; title2 'Residual plots for subsequent recoveries'; plot r_adj*p; data recov2; set recov; if p^=0 then lpred=log(p); else lpred=-10; proc plot data=recov2; title3 'Residual versus log(predicted)'; plot r_adj*lpred; run; data gof; set gof; p_value=1-CDF('CHISQ', D, df); proc print data=gof; title3 'Goodness of fit p-value'; run; quit;