************************************************************************* SAMPLE SUZE UPDATE COMPUTATIONS USING PARAMETRIC APPROACH C 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 ; /* Weibull Survival Function */ start WeiSurv(time,Wei1,Wei2) ; if time <= 0.0000001 then dum1=1; /* S(0)=1 */ else if time >= 1E99 then dum1=0 ; else dum1=exp(-Wei1 * (time**Wei2)) ; /* Median=[(1/Wei1) log 2]^(1/Wei2) */ return(dum1) ; finish ; /* ================================================================================================ */ /* ------------------ Jackup factor calculation ----------- */ /* For case: 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(7,7,0); CovMPhat[1:4,1:4]=CovPGhat ; CovMPhat[5,5]=pmissSE*pmissSE ; CovMPhat[6:7,6:7]=CovWehat ; 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 ; dpidweii=0 ; dpidweis=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=WeiSurv(dumtime[k,],Wei1est,Wei2est)-WeiSurv(dumtime[(k+1),],Wei1est,Wei2est) ; 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 ; if DPInd=0 then do; /* deals with case where Wei1hat=0 and Wei2hat=0 */ do k=(kl[i,]+1) to nt ; dum2=((1-pmisshat)##(kl[i,]-ml[i,]))*(pmisshat##(ml[i,]+k-kl[i,]-1)) ; dum4=0 ; dum5=0 ; dum6=0 ; if k>1 then do; dum4=(log(dumtime[k,])-WeiInt)/WeiSca ; dum5=exp(dum4) ; dum6=exp(-dum5) ; end ; dum7=0 ; dum8=0 ; dum9=0 ; if keps2)); *** 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 Weibull model to data, create information matrices and associated probabilities from model */ /* - 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 fit of Weibull model */ 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 */ DPInd=0 ; wei1est=0 ; wei2est=0 ; WeiInt=0 ; WeiSca=0 ; CovWehat=j(2,2,0) ; /* Initialize model estimates */ /* 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 LIFEREG if no dropout */ %if (&DPInd=1) %then %goto skipLife ; proc lifereg data=droprigc covout outest=dropestr noprint ; model ftime*fInd(1)= /d=weibull ; run ; /* Non-convergence will show up in OUTPUT */ proc iml workspace=512 symsize=256; load ; use dropestr ; read all into dum1 ; /* Get Weibull parameter estimates from right censored data */ WeiInt=dum1[1,2] ; /* Parameters in SAS parameterisation */ WeiSca=dum1[1,4] ; Wei1est=exp(-WeiInt/WeiSca) ; /* Parameter estimates in usual parameterisation */ Wei2est=1/WeiSca ; CovWehat=dum1[2:3,2]||dum1[2:3,4] ; /* Cov matrix for (SAS) parameters */ store _all_ module=_all_; quit ; /* Compute the prob associated with each observation pattern */ %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=WeiSurv(dumtime[k,],Wei1est,Wei2est)-WeiSurv(dumtime[(k+1),],Wei1est,Wei2est) ; 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=WeiSurv(dumtime[k,],Wei1est,Wei2est)-WeiSurv(dumtime[(k+1),],Wei1est,Wei2est) ; 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 ;