************************************************************ * * * One-way random effects model * * * * Artificial insemination of cows * * * * Outcome: Percentage of conceptions to * * services for successive samples * * * ************************************************************; /* Bulls example using method of moments (MM) approach */ OPTIONS NODATE NONUMBER; DATA bulls_summary; INPUT bull n mean var_mean; si_2 = n * var_mean; weight = n / si_2; weight_squared = weight**2; weight_mean = weight * mean; DATALINES; 1 5 41.2 35.14 2 2 64.5 30.25 3 7 56.3 18.99 4 5 39.6 101.06 5 7 67.1 38.64 6 9 53.2 27.72 ; /* Calculate some sums */ PROC SUMMARY DATA=bulls_summary SUM N; VAR weight weight_squared weight_mean; OUTPUT OUT=sums SUM=weight_sum weight_squared_sum weight_mean_sum N = k; RUN; /* Define some macro variables */ DATA _NULL_; SET sums; gd = weight_mean_sum / weight_sum; CALL SYMPUT ("gd", gd); CALL SYMPUT ("vhat_sum", weight_sum); CALL SYMPUT ("vhat_sum_sq", weight_squared_sum); CALL SYMPUT ("number",k); RUN; /* Compute Cochran's homogeneity statistic */ DATA cochran; SET bulls_summary; cochran_term = weight * (mean - &gd)**2; RUN; PROC SUMMARY DATA=cochran SUM; VAR cochran_term; OUTPUT OUT=cochran_sum SUM=cochran_sum; RUN; /* DerSimonian-Laird estimator */ DATA _NULL_; SET cochran_sum(KEEP=cochran_sum); tau_dsl = (cochran_sum - (&number - 1)) / (&vhat_sum - &vhat_sum_sq / &vhat_sum); tau_dsl = max(0, tau_dsl); CALL SYMPUT ("DSL",tau_dsl); RUN; /* Random effects meta-analysis */ DATA random; SET bulls_summary; weight_ran = 1 / (si_2 / n + &DSL); weight_ran_mean = weight_ran * mean; RUN; PROC SUMMARY DATA=random SUM; VAR weight_ran weight_ran_mean; OUTPUT OUT=random_sums SUM=weight_ran_sum weight_ran_mean_sum; RUN; DATA _NULL_; SET random_sums; est_random = weight_ran_mean_sum / weight_ran_sum; CALL SYMPUT ("est_ran",est_random); CALL SYMPUT ("what_sum",weight_ran_sum); RUN; /* Calculating Hartung's variance estimator */ DATA hartung; SET bulls_summary; weight_ran = 1 / (si_2 / n + &DSL); q_term = weight_ran *(mean - &est_ran)**2; RUN; PROC SUMMARY DATA=hartung SUM; VAR q_term; OUTPUT OUT=hartung_sum SUM=q_sum; RUN; DATA _NULL_; SET hartung_sum; q_sum = q_sum / (&what_sum * (&number - 1)); CALL SYMPUT("q_sum",q_sum); RUN; /* Random effects meta-analysis */ DATA random_MA; LENGTH method $40.; method = "FE Meta-Analysis"; estimate = &gd; lower_CI = &gd - sqrt(1/&vhat_sum) * QUANTILE("NORMAL",0.975,0,1); upper_CI = &gd + sqrt(1/&vhat_sum) * QUANTILE("NORMAL",0.975,0,1); OUTPUT; method = "RE Meta-Analysis (DerSimonian-Laird)"; estimate = &est_ran; lower_CI = &est_ran - sqrt(1/&what_sum) * QUANTILE("NORMAL",0.975,0,1); upper_CI = &est_ran + sqrt(1/&what_sum) * QUANTILE("NORMAL",0.975,0,1); OUTPUT; method = "RE Meta-Analysis (Hartung-Knapp)"; estimate = &est_ran; lower_CI = &est_ran - sqrt(&q_sum) * QUANTILE("T",0.975,&number - 1); upper_CI = &est_ran + sqrt(&q_sum) * QUANTILE("T",0.975,&number - 1); OUTPUT; RUN; PROC REPORT DATA=random_MA HEADLKNE HEADSKIP NOWINDOWS SPLIT="*"; COLUMNS method estimate lower_CI upper_CI; DEFINE method / DISPLAY "Method"; DEFINE estimate / DISPLAY "Meta-Analysis*Estimate" FORMAT = 15.4; DEFINE lower_CI / DISPLAY "95% CI Lower Bound" FORMAT = 15.4; DEFINE upper_CI / DISPLAY "95% CI Upper Bound" FORMAT = 15.4; TITLE "Artificial Insemination Example"; RUN;