This script requires sas/stat 13.1 or newer to perform repeated measures using proc glmpower. Code below checks to verify version.
PROC PRODUCT_STATUS;
run;
Location of data files
libname files "/mnt/hgfs/myfolders";
The following macro was written and employed to perform the sums of squares calculation required to estimate the residual variance data from the pilot study.
%macro SS (A,B);
%let output_title = &B;
PROC iml;
use &A;
read all var {head body tail} into x;
read all var {head} into head;
read all var {body} into body;
read all var {tail} into tail;
meanhead = mean(head);
meanbody = mean(body);
meantail = mean(tail);
rowMeans = x[ ,:];
colMeans = x[:, ];
grandMean= x[:];
*to calculate SSw (within-group SS);
part1=ssq(head - meanhead);
part2=ssq(body-meanbody);
part3=ssq(tail-meantail);
SSw=(part1+part2+part3);
*to calculate SSs (within subject SS);
part4=ssq(rowmeans[1]-grandmean);
part5=ssq(rowmeans[2]-grandmean);
part6=ssq(rowmeans[3]-grandmean);
SSs=3*(part4+part5+part6);
*to calculate SSe (error term);
SSe=SSw-SSs;
*to get SSe variance, where 2 is for degrees of freedom;
var=Sse/2;
*to get SD;
SD=sqrt(var);
*to get SScondition;
SScondition=3*(ssq(meanhead-grandmean)+ssq(meanbody-grandmean)+ssq(meantail-grandmean));
*to get SStotal (total SS);
SStotal=SSw+SScondition;
print SStotal SScondition SSw SSs SSe var SD;
create SS_&B var {"SStotal" "SScondition" "SSw" "SSs" "SSe" "var" "SD"};
append;
close SS_&B;
quit;
title "&output_title";
%mend ss;
quit;
To perform the power calculations, an estimate of the residual variance (i.e. error variability, where error plus subject variability equals total within-group variability) for each of the three markers was calculated using pilot data from nPOD donors who were NOT used in the main study. The process includes data 1) import, 2) transformation, and 3) analysis.
For CD45, SD was found to be 1.98%
*import;
DATA cd45;
set files.cd45;
run;
DATA cd45_pilot;
set cd45;
if case^=6222 and case^=6092 and case^=6107 then delete; *only need pilot data;
if case=6222 then match_group=1;
if case=6107 then match_group=2;
if case=6092 then match_group=3;
run;
*transformation;
PROC sort data=cd45_pilot;
by case match_group hospbin;
run;
PROC transpose data=cd45_pilot out=cd45_pilot1 (rename=( a=Head b=Body c=Tail));
by case match_group hospbin;
id sample_type2;
var percent_cd45;
run;
*analysis;
%ss(cd45_pilot1, cd45_ss);
For CD68, SD was found to be 0.55%
*import;
DATA cd68;
set files.cd68;
run;
DATA cd68_pilot;
set cd68;
if case^=6222 and case^=6092 and case^=6107 then delete; *only need pilot data;
if case=6222 then match_group=1;
if case=6107 then match_group=2;
if case=6092 then match_group=3;
run;
*transformation;
PROC sort data=cd45_pilot;
by case match_group hospbin;
run;
PROC transpose data=cd68_pilot out=cd68_pilot1 (rename=( a=Head b=Body c=Tail));
by case match_group hospbin;
id sample_type2;
var percent_cd68;
run;
*analysis;
%ss(cd68_pilot1, cd68_ss);
For Ki67, SD was found to be 0.74%
*import;
DATA ki67;
set files.ki67;
run;
DATA ki67_pilot;
set ki67;
if case^=6222 and case^=6092 and case^=6107 then delete; *only need pilot data;
if case=6222 then match_group=1;
if case=6107 then match_group=2;
if case=6092 then match_group=3;
run;
*transformation;
PROC sort data=ki67_pilot;
by case match_group hospbin;
run;
PROC transpose data=ki67_pilot out=ki67_pilot1 (rename=( a=Head b=Body c=Tail));
by case match_group hospbin;
id sample_type2;
var percent_ki67;
run;
*analysis;
%ss(ki67_pilot1, ki67_ss);
Short (<3 days) duration hospitalization group: To obtain theoretical values: 1) Cell number fractions for islet (0.483), acinar (0.253), duct (0.227), and other cell types (0.038; includes leukocytes, neuronal elements, dead, stromal, and endothelial cells) found in an islet prep were obtained from Pisania et al (2010, Laboratory Investigation); it was assumed that leukocytes contributed 1/5th of the other cell types (i.e. 0.0076), 2) The total cell number fraction was adjusted from 1.001 to 0.774 (i.e. 1.001 minus ductal fraction, since this cell type was not analyzed in the experiment), 3) The individual cell number fractions, minus duct, for leukocytes (0.0076/0.774=0.0098) and islets (0.483/0.774=0.624) was calculated, 4) The number of leukocytes/islet was calculated (0.098/0.624=0.015) in islet preps and compared to 0.1 $\pm$ 0.1 for whole tissue (Bogdani et al, 2014, Diabetes). Therefore, the initial estimate of 0.98% (step 3) was conservatively doubled to 1.96% to account for the possibility of an increased % of leukocytes in whole tissue.
Mid (>=3 to < 6 days) and long (>=6 days) duration hospitalization groups: See main text. In brief, the magnitude of the differences in number of cells counted by In'tVeld and colleagues (2010, Diabetes) were used (Table 5, <=25 year old group) to adjust up the estimates from the short duration group.
Differences in leukocyte % based on pancreas region: The authors are not aware of any data suggesting differences based on region of pancreas, therefore it was assumed that % remained constant in the head, body, and tail, i.e. 1.96%.
Step 1: Creation of exemplary dataset (mean values)
DATA cd45;
input hospbin $ head body tail alloc;
datalines;
ICU_short 1.96 1.96 1.96 1
ICU_mid 1.96 1.96 1.96 1
ICU_long 5.88 5.88 5.88 1
;
run;
Step 2: Additional assumptions and analysis
PROC glmpower data = cd45;
class hospbin;
model head body tail = hospbin;
repeated panc_region;
weight alloc;
power
mtest=uncorr gg
/*available options include uncorr or GG*/
/*uncorr used as analysis plan is to report unadjusted univariate p values for within and between subject effects, including interaction term*/
/*If sphericity is violated, GG will be used for univariate within-subjects effects*/
stddev = 1.9800816
/*data taken from pilot study*/
matrix ("mycorrmat")= lear(0.99, 0, 3, 1 2 3)
/*base corr, corr-decay, number of levels, level-values (defaul is used if left blank)*/
/*base corr and corr-decay specified as constant correlation of 0.99 with no decay, i.e. clustered correlation or compound symmetry*/
/*0.99 calculated as the mean of head body, body tail, and head tail correlations of panc weight in control donors without missing data*/
corrmat = "mycorrmat"
alpha = 0.05
ntotal =39
power = .;
/*option to plot x=power min = .6 max =1*/
run;
Step 1: Creation of exemplary dataset (mean values)
DATA cd68;
input hospbin $ head body tail alloc;
datalines;
ICU_short 1.96 1.96 1.96 1
ICU_mid 2.94 2.94 2.94 1
ICU_long 5.88 5.88 5.88 1
;
run;
Step 2: Additional assumptions and analysis
PROC glmpower data = cd68;
class hospbin;
model head body tail = hospbin;
repeated panc_region;
weight alloc;
power
mtest=uncorr gg
/*available options include uncorr or GG*/
/*uncorr used as analysis plan is to report unadjusted univariate p values for within and between subject effects, including interaction term*/
/*If sphericity is violated, GG will be used for univariate within-subjects effects*/
stddev = 0.5488278
/*data taken from pilot study*/
matrix ("mycorrmat")= lear(0.99, 0, 3, 1 2 3)
/*base corr, corr-decay, number of levels, level-values (defaul is used if left blank)*/
/*base corr and corr-decay specified as constant correlation of 0.99 with no decay, i.e. clustered correlation or compound symmetry*/
/*0.99 calculated as the mean of head body, body tail, and head tail correlations of panc weight in control donors without missing data*/
corrmat = "mycorrmat"
alpha = 0.05
ntotal =39
power = .;
/*option to plot x=power min = .6 max =1*/
run;
Step 1: Creation of exemplary dataset (See main text for derivation of initial values).
DATA ki67;
input hospbin $ head body tail alloc;
datalines;
ICU_short 0.1 0.1 0.1 1
ICU_mid 0.6 0.6 0.6 1
ICU_long 1.59 1.59 1.59 1
;
run;
Step 2: Additional assumptions and analysis
PROC glmpower data = ki67;
class hospbin;
model head body tail = hospbin;
repeated panc_region;
weight alloc;
power
mtest=uncorr gg
/*available options include uncorr or GG*/
/*uncorr used as analysis plan is to report unadjusted univariate p values for within and between subject effects, including interaction term*/
/*If sphericity is violated, GG will be used for univariate within-subjects effects*/
stddev = 0.7407297
/*data taken from pilot study*/
matrix ("mycorrmat")= lear(0.99, 0, 3, 1 2 3)
/*base corr, corr-decay, number of levels, level-values (defaul is used if left blank)*/
/*base corr and corr-decay specified as constant correlation of 0.99 with no decay, i.e. clustered correlation or compound symmetry*/
/*0.99 calculated as the mean of head body, body tail, and head tail correlations of panc weight in control donors without missing data*/
corrmat = "mycorrmat"
alpha = 0.05
ntotal =39
power = .;
/*option to plot x=power min = .6 max =1*/
run;