************************************************************************* SAMPLE SIZE UPDATE COMPUTATIONS USING NONPARAMETRIC APPROACH B Zucker, D. M. and Denne, J. (2002). Sample size redetermination for repeated measures studies. Biometrics, to appear. ************************************************************************* ************************************************************************* Copyright, 22 May 2002, D Zucker and J Denne. All rights reserved. NOTE TO USERS: This code is released on an "AS IS" basis. THERE IS NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED. Users are warned that the developers are not liable for any errors or damages that occur as a result of using this code. If you find these programs of help in your research we ask you please to acknowledge their source in reports and publications. Permission to use, copy, modify, and distribute this software and its documentation for any purpose, without fee, and without a written agreement, is hereby granted, provided that the above copyright notice and this NOTE TO USERS appears in all copies. *************************************************************************; options ls=80 nodate ps=55 nocenter; /*--------------------------------------------------------------------------------------------------- */ /*--------------------------------------------------------------------------------------------------- */ /*----------------------- SUBROUTINES -------------------------------------------------------------- */ proc iml workspace=512 symsize=256; /* Short function used in creating matpat */ /* the matrix containing all possible missingness patterns */ /* Method: move along line, Find a zero, change to a one, then stop and return */ /* Find a 1, change to a zero and continue along line */ start mpt(x,m) ; do j=1 to m ; if(x[,j] > 0.5 ) then goto t3 ; x[,j]=x[,j]+1 ; return(x) ; t3:x[,j]=0 ; end ; return(x) ; finish ; /* ================================================================================================ */ /* ------------------ Jackup factor calculation ----------- */ /* For case: KM non-parametric estimate */ /* FOR INPUTS SEE EXAMPLE CODE BELOW */ start jack; /* Form asymptotic cov matrix for covariance/missingness parameters */ /* Have to work in SAS parameterisation for Weibull parameters */ use covdum ; read all into CovPGhat ; CovPGhat=CovPGhat[1:4,2:5] ; CovMPhat=j(5+nt-1,5+nt-1,0); CovMPhat[1:4,1:4]=CovPGhat ; CovMPhat[5,5]=pmissSE*pmissSE ; CovKM=j(nt,nt,0); do r=1 to nt; var=surv[r,3]**2; CovKM[r,r]=var; do s=r+1 to nt; CovKM[r,s] = var; CovKM[s,r] = var; end; end; CovMPhat[(6:6+nt-2),(6:6+nt-2)] = CovKM[(2:nt),(2:nt)] ; /* cov matrix of post-baseline KM surv estimates */ dumtime=time[(1:nt),]//1E99 ; do i=1 to totpat ; dumt=time[loc(matpat[i,]=1),] ; zij=j(nrow(dumt),1,1)||dumt ; eij=I(nrow(dumt))*errvhat ; /* R- matrix assumed to be Identity */ invsig=inv(zij*Ghat*t(zij)+eij) ; /* create inverse of Sigma_ij */ xdum=j(nrow(dumt),1,1)||j(nrow(dumt),1,0)||dumt||j(nrow(dumt),1,0); /* Xmatrix for subject from trt group 0 */ Ui=t(xdum)*invsig*xdum ; dUi1=-t(xdum)*invsig*zij*{1 0, 0 0}*t(zij)*invsig*xdum ; dUi2=-t(xdum)*invsig*zij*{0 1, 1 0}*t(zij)*invsig*xdum ; dUi3=-t(xdum)*invsig*zij*{0 0, 0 1}*t(zij)*invsig*xdum ; dUi4=-t(xdum)*invsig*invsig*xdum ; xdum=j(nrow(dumt),1,0)||j(nrow(dumt),1,1)||j(nrow(dumt),1,0)||dumt; /* Xmatrix for subject from trt group 1 */ Ui=(Ui+t(xdum)*invsig*xdum)/2 ; dUi1=probpat[i,]*(dUi1-t(xdum)*invsig*zij*{1 0, 0 0}*t(zij)*invsig*xdum)/2 ; dUi2=probpat[i,]*(dUi2-t(xdum)*invsig*zij*{0 1, 1 0}*t(zij)*invsig*xdum)/2 ; dUi3=probpat[i,]*(dUi3-t(xdum)*invsig*zij*{0 0, 0 1}*t(zij)*invsig*xdum)/2 ; dUi4=probpat[i,]*(dUi4-t(xdum)*invsig*invsig*xdum)/2 ; dpidpmis=0 ; if (pmisshat >= 0.00000001 | pmisshat <= 0.99999999) then do ; /* deals with case where pmisshat=0 or 1 */ do k=(kl[i,]+1) to nt ; dum1=Surv[k,2]-Surv[(k+1),2] ; ; dum2=((1-pmisshat)##(kl[i,]-ml[i,]))*(pmisshat##(ml[i,]+k-kl[i,]-1)) ; dum3=(((ml[i,]+k-kl[i,]-1)/pmisshat)-((kl[i,]-ml[i,])/(1-pmisshat))) ; dpidpmis=dpidpmis+dum1*dum2*dum3 ; end ; end ; dpids=j((nt-1),1,0) ; /* derivative of pi wrt postbaseline KM surv estimates */ if DPInd=0 then do; /* does not do if no drop-out */ do k=1 to (nt-1) ; /* derivative of pi wrt postbaseline KM surv estimates */ if k>=(kl[i,]+1) then dpids[k,1]=-((1-pmisshat)##(kl[i,]-ml[i,]))*(pmisshat##(ml[i,]+k-kl[i,]-1)) ; if (k<=(nt-1) & k >=kl[i,]) then dpids[k,1]=dpids[k,1]+((1-pmisshat)##(kl[i,]-ml[i,]))*(pmisshat##(ml[i,]+k-kl[i,])) ; end ; end ; dUi5=dpidpmis*Ui ; dUi6=dpids[1,1]*Ui ; dUi7=dpids[2,1]*Ui ; dUi8=dpids[3,1]*Ui ; dUi9=dpids[4,1]*Ui ; if i=1 then dH1=dUi1 ; else dH1=dH1+dUi1 ; if i=1 then dH2=dUi2 ; else dH2=dH2+dUi2 ; if i=1 then dH3=dUi3 ; else dH3=dH3+dUi3 ; if i=1 then dH4=dUi4 ; else dH4=dH4+dUi4 ; if i=1 then dH5=dUi5 ; else dH5=dH5+dUi5 ; if i=1 then dH6=dUi6 ; else dH6=dH6+dUi6 ; if i=1 then dH7=dUi7 ; else dH7=dH7+dUi7 ; if i=1 then dH8=dUi8 ; else dH8=dH8+dUi8 ; if i=1 then dH9=dUi9 ; else dH9=dH9+dUi9 ; end ; * print dUi5 dUi6 dUi7 dUi8 dUi9 ; delh=j(9,10,999) ; /* Form Del h */ ind=0 ; do sk=1 to 4 ; /* Work down columns of dH from lead diagonal */ do sj=sk to 4 ; ind=ind+1 ; delh[,ind]=dH1[sj,sk]//dH2[sj,sk]//dH3[sj,sk]//dH4[sj,sk]//dH5[sj,sk]//dH6[sj,sk]//dH7[sj,sk]//dH8[sj,sk]//dH9[sj,sk] ; end ; end ; covhhat=t(delh)*CovMPhat*delh ; /* Covariance matrix for vech(h) */ /* Now calculate H-hat and alpha-hat */ ahat=inv(Hhat)*c ; dum1=ahat#ahat ; dum2=2*ahat*t(ahat)-diag(dum1) ; alphhat=dum2[1:4,1]//dum2[2:4,2]//dum2[3:4,3]//dum2[4,4] ; /* Now to compute the estimated variance of tau^2 */ vartau2=t(alphhat)*covhhat*alphhat ; /* Now to compute the jack-up factor */ chidf=2*tau2*tau2/vartau2 ; if chidf < 2 then chidf=2 ; jackmult=jackchi(chidf,alpha,beta) ; finish jack; /* ================================================================================================ */ /* Jack-up calculation: integral over chi-sq */ start jackchi(nu,alpha,beta) ; *** nu= DEGREES OF FREEDOM ***; *** SETUP ***; sq2pi = sqrt(2*arcos(-1)); pow = 1 - beta; za = probit(1-alpha); zb = probit(pow); zfac = za + zb; maxiter = 20; n_part = 500; eps1 = 0.000001; eps2 = 0.0001; h = (1-eps1) / n_part; m = 1; om = 1.0; diff = 99; *** NEWTON ITERATION LOOP ***; do while ((m<=maxiter) & (diff>eps2)); *** NUMERICAL INTEGRATION LOOP ***; power = 0; pder = 0; do i = 0 to n_part; wt = 2 * (1+mod(i,2)); if ((i=0) | (i=n_part)) then wt=1; u1 = i*h; v1 = cinv(u1,nu)/nu; s1 = sqrt(v1); arg = zfac*s1*om-za; power = power + wt*probnorm(arg); pder = pder + wt*zfac*s1*exp(-0.5*(arg**2)); end; power = (h/3) * power; pder = (h/3) * pder / sq2pi; omold = om; om = omold - (power-pow)/pder; diff = abs(om-omold); end; om2 = om**2; return(om2) ; finish jackchi ; store _all_ module=_all_; quit ; /* ================================================================================================ */ /* ================================================================================================ */ /* EXAMPLE INPUT TO subroutine JACK to calculate Inflation Factor */ /* Internal Pilot data : include all subjects who have data at the time of the internal pilot */ /* Example data */ data data1 ; input SUBJECT TREAT TIME RESPONSE; cards ; 1 0 0.00 -1.52650 1 0 0.25 -1.04849 1 0 0.50 -0.47514 1 0 0.75 -0.47329 1 0 1.00 0.99631 2 1 0.00 -1.74538 2 1 0.25 -1.09840 2 1 0.50 -1.93634 2 1 0.75 -1.96460 2 1 1.00 -1.70901 3 0 0.00 0.90776 3 0 0.25 0.18904 3 0 0.50 1.98844 3 0 0.75 1.88236 3 0 1.00 3.86319 4 1 0.00 -0.72858 4 1 0.25 -0.50691 5 0 0.00 0.96746 5 0 0.25 0.78068 5 0 0.50 0.70842 5 0 0.75 0.59748 6 1 0.00 -0.75216 7 0 0.00 -1.72473 8 1 0.00 0.18941 8 1 0.25 -0.67676 8 1 0.50 0.99921 8 1 0.75 0.51826 8 1 1.00 0.85929 9 0 0.00 -0.54776 9 0 0.25 -0.52594 9 0 0.50 -0.06257 9 0 0.75 -0.00790 9 0 1.00 0.55285 10 1 0.00 -1.21507 10 1 0.25 -1.17854 10 1 0.50 -1.88146 10 1 0.75 -1.87816 10 1 1.00 -2.43429 11 0 0.00 -0.27668 11 0 0.25 -1.24699 11 0 0.50 -2.42856 12 1 0.00 0.32682 12 1 0.25 0.04770 13 0 0.00 1.62511 13 0 0.25 1.52302 13 0 0.50 1.40110 13 0 0.75 3.29042 14 1 0.00 1.48042 14 1 0.25 2.34662 14 1 0.50 3.46417 14 1 0.75 4.19304 15 0 0.00 -0.35491 15 0 0.25 -0.63751 15 0 0.50 -0.47484 15 0 0.75 -0.60756 16 1 0.00 0.89835 16 1 0.25 0.83360 16 1 0.50 0.63680 16 1 0.75 0.58950 17 0 0.00 1.95192 17 0 0.25 3.86006 17 0 0.50 3.98230 18 1 0.00 -1.04555 18 1 0.25 -1.30094 18 1 0.50 -1.89748 19 0 0.00 -0.17763 19 0 0.25 0.57057 19 0 0.50 0.52868 20 1 0.00 -0.26452 20 1 0.25 0.65028 20 1 0.50 -0.12757 21 0 0.00 -1.23502 22 1 0.00 -1.40183 22 1 0.25 -0.68828 23 0 0.00 -1.02030 24 1 0.00 -1.37771 24 1 0.25 -0.51025 25 0 0.00 0.97804 26 1 0.00 -1.98672 27 0 0.00 0.31216 28 1 0.00 -0.19089 29 0 0.00 0.00418 30 1 0.00 0.17965 ; run; /* Get estimates of covariance parameters */ proc mixed asycov noclprint noitprint data=data1; class treat ; model response=treat time treat*time /s ddfm=satterth ; random int time /sub=subject type=un ; make 'Covparms' out=gdum noprint ; make 'Asycov' out=covdum noprint ; make 'Fitting' out=garb1 noprint ; make 'Tests' out=garb2 noprint ; make 'SolutionF' out=Soldum1 noprint ; run ; proc iml workspace=512 symsize=256; load ; alpha=0.025 ; beta=0.1 ; /* One-sided Type I and Power */ delta=1.178 ; /* Difference to be detected with power 1-beta */ nt=5 ; /* maximum number of (equally-spaced) visits */ /* written for nt=5 - code for estimating survival will need minor alteration if nt is changed */ time=t((1:(2*nt-4))-1)/(nt-1) ; /* Planned Visit times: nt originally planned, times <=1. Inc baseline */ use data1 ; read all var{subject} into dum1 ; nsub=dum1[nrow(dum1),] ; /* Read in number of subjects (nsub); */ /* Put estimated Gmatrix and error variance into an IML matrix and variable */ use gdum ; read all into dum1 ; Ghat=j(2,2,1) ; Ghat[1,1]=dum1[1,] ; Ghat[1,2]=dum1[2,] ; Ghat[2,1]=dum1[2,] ; Ghat[2,2]=dum1[3,] ; errvhat=dum1[4,] ; store _all_ module=_all_; quit ; /* Internal Pilot data : include all subjects who have data at the time of the internal pilot */ /* Example data */ /* FTIME = drop-out time or censoring time */ /* FIND = indicator, 0=drop-out 1=censored */ data droprigc ; input SUBJECT FTIME FIND ; cards ; 1 1.00000 1 2 1.00000 1 3 1.00000 1 4 0.30907 0 5 0.82787 0 6 0.00336 0 7 0.13334 0 8 1.00000 1 9 1.00000 1 10 1.00000 1 11 0.52606 0 12 0.25184 0 13 0.86842 1 14 0.86842 1 15 0.76316 1 16 0.76316 1 17 0.65789 1 18 0.65789 1 19 0.55263 1 20 0.55263 1 21 0.24229 0 22 0.44737 1 23 0.01409 0 24 0.34211 1 25 0.23684 1 26 0.23684 1 27 0.13158 1 28 0.13158 1 29 0.02632 1 30 0.02632 1 run ; %macro fitest ; /* Macro to fit KM curve to data, create information matrices and associated probabilities */ /* - note we assume that baseline cannot be missing */ /* == NOTE: code for calculating estimates for intermittent missingness (pmiss and pmissSE) not provided */ /* == entered manually here */ /* Get estimates of survival from Kaplan Meier curve */ proc iml workspace=512 symsize=256; load; /* Intermittent missingness - note we assume that baseline cannot be missing */ nmiss=0 ; /* Hand calculation of number of intermittently missed visits prior to drop-out/censoring in internal pilot */ nvis=53 ; /* Hand calculation of maximum possible number of visits prior to drop-out/censoring in internal pilot */ pmisshat=nmiss/nvis ; /* proportion of visits missed by patients prior to drop-out */ pmissSE=(pmisshat*(1-pmisshat)/nvis)##(0.5) ; /* SE for proportion of visits missed by patients prior to drop-out */ surv={0 1 0,0.25 1 0,0.5 1 0,0.75 1 0,1.0 1 0, 1E99 0 0}; /* set up initial value incase no missing data */ DPInd=0 ; /* Deal with case where no observed dropouts */ use droprigc ; read all var{subject ftime fInd} into xydum ; DPInd=all(xydum[,3]) ; /* Indicator variable, =1 if no observed events for right censored data */ create Droppar from DPInd[ colname={Ind}]; append from DPInd ; /* Make macro variable indicator skipping LIFEREG */ store _all_ module=_all_; quit ; Data Droppar; set Droppar ; call symput('DPInd',Ind) ; run ; /* Make macro variable indicator skipping LIFETEST if no dropout */ %if (&DPInd=1) %then %goto skipLife ; data droprigc ; set droprigc ; /* change all censored at time=1 to obs fail at time=5 */ ftime1=ftime ; fInd1=fInd ; /* forces SAS to give estimate at time=1 */ if (ftime=1 & fInd=1) then do ; ftime1=5 ; fInd1=0 ; end ; run ; proc lifetest data=droprigc outs=surv timelist=(0.25,0.5,0.75,1) reduceout noprint ; time ftime1*fInd1(1) ; run ; data surv ; set surv ; sdkm=(survival-SDF_LCL)/1.959964 ; keep _time_ survival sdkm ; run ; proc iml workspace=512 symsize=256; load ; use surv ; read all into surv ; surv={0 1 0}//surv//{1E99 0 0} ; do k=2 to nt ; if surv[k,2]=. then surv[k,2]=surv[k-1,2] ; if surv[k,3]=. then surv[k,3]=surv[k-1,3] ; end ; store _all_ module=_all_; quit ; /* Compute information matrix for each visit pattern and associated probability */ %skipLife: proc iml workspace=512 symsize=256; load ; if pmisshat=0 then pmisshat=0.000000001 ; if pmisshat=1 then pmisshat=0.999999999 ; totpat=2**(nt-1) ; matpat=j(totpat,(nt-1),999) ; /* Create matrix of possible patterns of missingness */ probpat=j(totpat,1,0) ; /* Assoc probbility of missingness patern */ kl=j(totpat,1,0) ; ml=j(totpat,1,0) ; matpat[1,]=j(1,(nt-1),0) ; dumtime=time[(1:nt),]//1E99 ; do k=(kl[1,]+1) to nt ; dum1=Surv[k,2]-Surv[(k+1),2] ; probpat[1,]=probpat[1,]+dum1*((1-pmisshat)##(kl[1,]-ml[1,])*pmisshat##(ml[1,]+k-kl[1,]-1)) ; end ; do i=2 to totpat ; dum1=matpat[i-1,] ; matpat[i,]=mpt(dum1,(nt-1)) ; kl[i,]=1+int(log(i-0.0001)/log(2)) ; /* kl is position of last 1 in row of matpat before adding col of 1's */ ml[i,]=kl[i,]-sum(matpat[i,]) ; /* ml is number of zero's before last 1 in each row */ do k=(kl[i,]+1) to nt ; dum1=Surv[k,2]-Surv[(k+1),2] ; probpat[i,]=probpat[i,]+dum1*((1-pmisshat)##(kl[i,]-ml[i,])*pmisshat##(ml[i,]+k-kl[i,]-1)) ; end ; end ; matpat=j(totpat,1,1)||matpat ; /* baseline never missing */ dum3=sum(probpat) ; if dum3 <0.9999999 | dum3 > 1.000001 then print dum3 ; /* check that probabilities sum to 1 */ /* Now to use the patterns and assoc probs to compute tau2 */ do i=1 to totpat ; dumt=time[loc(matpat[i,]=1),] ; zij=j(nrow(dumt),1,1)||dumt ; eij=I(nrow(dumt))*errvhat ; /* R- matrix assumed to be Identity */ invsigij=inv(zij*Ghat*t(zij)+eij) ; /* create inverse of Sigma_ij */ xdum0=j(nrow(dumt),1,1)||j(nrow(dumt),1,0)||dumt||j(nrow(dumt),1,0); /* Xmatrix for subject from trt group 0 */ xdum1=j(nrow(dumt),1,0)||j(nrow(dumt),1,1)||j(nrow(dumt),1,0)||dumt; /* Xmatrix for subject from trt group 1 */ if i=1 then Hhat=(t(xdum0)*invsigij*xdum0+t(xdum1)*invsigij*xdum1)*probpat[i,]/2 ; /* Compute Information matrix */ else Hhat=Hhat+(t(xdum0)*invsigij*xdum0+t(xdum1)*invsigij*xdum1)*probpat[i,]/2 ; end ; if det(Hhat)=0 then print Hhat ; /* check for non-signularity in Hhat */ c={0,0,-1,1} ; /* Define constrast */ tau2=t(c)*inv(Hhat)*c ; store _all_ module=_all_; quit ; %mend ; %fitest ; /* calculate jackmult, the inflation factor and hence the sample size */ proc iml workspace=512 symsize=256; load ; jackmult=1 ; run jack ; /* FIX TARGET SAMPLE SIZE PER ARM */ mt1=jackmult*0.5*tau2*(probit(1-alpha)+probit(1-beta))**2/(delta**2) ; print alpha beta delta chidf jackmult mt1 ; store _all_ module=_all_; quit ;