************************************************************************* INFLATION FACTOR COMPUTATIONS FOR NONPARAMETRIC APPROACH A 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; proc iml; load ; /* ================================================================================================ */ /* ------------------ Jackup factor calculation for Non-parametric SSR ----------- */ /* INPUTS: */ /* Data-set data1: structure - variables in following order */ /* SUBJECT = subject number */ /* TREAT = treatment (0 or 1) */ /* TIME = time of observation or visit number */ /* */ /* errvhat = variance of error term (output from PROC MIXED) */ /* Ghat = 2x2 covariance matrix for random intercepts/slopes (output from PROC MIXED) */ /* CovPGhat= asymptotic cov matrix of covariance parameters */ /* nsub = number of subjects in data1 */ /* tau2 = nsub * variance of contrast of interest from model fit to data1 */ /* alpha = One-sided Type I */ /* beta = Type II error rate */ /* */ /* */ start jack; do i=1 to nsub; use data1 ; read all var{subject treat time} where (subject=i) into x ; /* Read in data for each subject sequentially */ Rmat=I(nrow(x)) ; /* R- matrix taken to be Identity, but could be estimated */ dum1= j(nrow(x),1,1) ; zij=dum1||x[,3] ; eij=Rmat*errvhat ; sigmaij=zij*Ghat*t(zij)+eij ; /* create covariance matrix for subject i, Sigma_ij */ /* create (expanded) X */ dum2=x[,2] ; dum1=j(nrow(x),1,1)-x[,2] ; /* Indicators for treatment */ dum4=dum2#x[,3] ; dum3=dum1#x[,3] ; /* Coeff's for slope */ xdum=dum1||dum2||dum3||dum4 ; invsig=inv(sigmaij) ; 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 ; if i=1 then Umat=(i*j(4,1,1))||Ui ; else Umat=Umat//((i*j(4,1,1))||Ui) ; /* Store Ui's in Umat */ if i=1 then dUmat1=dUi1 ; else dUmat1=dUmat1+dUi1 ; /* Compute derivative matrix */ if i=1 then dUmat2=dUi2 ; else dUmat2=dUmat2+dUi2 ; /* Compute derivative matrix */ if i=1 then dUmat3=dUi3 ; else dUmat3=dUmat3+dUi3 ; /* Compute derivative matrix */ if i=1 then dUmat4=dUi4 ; else dUmat4=dUmat4+dUi4 ; /* Compute derivative matrix */ if i=1 then Inf=Ui ; else Inf=Inf+Ui ; /* Compute Information matrix */ end ; close data1 ; /* Now calculate H-hat and alpha-hat */ c={0,0,-1,1} ; Hhat=Inf/(nsub) ; /* H-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] ; /* Get asymptotic cov matrix of covariance parameters from PROC MIXED fit on internal pilot data */ use covdum ; read all into CovPGhat ; CovPGhat=CovPGhat[1:4,2:5] ; /* Now to calculate covariances of the elements of H with (sj,sk) and (sl,sm) being the (j,k)th and (l,m)th elements */ covhhat=j(10,10,0) ; rowind=(1:10) || (2:10) || (3:10) || (4:10) || (5:10) || (6:10) || (7:10) || (8:10) || (9:10) || 10 ; colind={ [10] 1 [9] 2 [8] 3 [7] 4 [6] 5 [5] 6 [4]7 [3] 8 [2] 9 [1] 10 } ; ind=0 ; do sk=1 to 4 ; /* Work down columns from lead diagonal */ do sj=sk to 4 ; do sm=sk to 4 ; dum1=sj ; if(sm > sk) then dum1=sm ; do sl=dum1 to 4 ; ind=ind+1 ; covtot=0 ; do i=1 to nsub; /* Loop over subjects */ Ui=Umat[loc(Umat[,1]=i),] ; Ui=Ui[1:4,2:5] ; /* remove subject label */ covtot=covtot+(Ui[sj,sk]-Hhat[sj,sk])*(Ui[sl,sm]-Hhat[sl,sm]) ; end ; dUjk=(dUmat1[sj,sk]//dUmat2[sj,sk]//dUmat3[sj,sk]//dUmat4[sj,sk])/(nsub) ; dUlm=(dUmat1[sl,sm]//dUmat2[sl,sm]//dUmat3[sl,sm]//dUmat4[sl,sm])/(nsub) ; covhhat[rowind[ind],colind[ind]]=((covtot/(nsub))+t(dUjk)*(nsub)*CovPGhat*dUlm)/(nsub) ; covhhat[colind[ind],rowind[ind]]=covhhat[rowind[ind],colind[ind]] ; end ; end ; end ; end ; /* 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 ; 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 completed, or would have completed had they not discontinued */ data data1 ; input SUBJECT TREAT TIME RESPONSE ; cards ; 1 0 0.00 -0.51968 1 0 0.25 -0.90984 1 0 0.50 -0.83201 1 0 0.75 -0.92290 1 0 1.00 -0.89301 2 1 0.00 -0.03281 3 0 0.00 -2.08366 3 0 0.25 -2.65272 3 0 0.50 -1.83465 3 0 0.75 -2.34985 3 0 1.00 -1.93915 4 1 0.00 -0.08825 4 1 0.25 -0.53483 4 1 0.50 -0.37673 4 1 0.75 -1.84424 5 0 0.00 0.12560 5 0 0.25 0.82496 5 0 0.50 0.78801 5 0 0.75 1.05155 5 0 1.00 2.13416 6 1 0.00 0.18603 7 0 0.00 0.42667 7 0 0.25 -0.26634 7 0 0.50 -1.05864 7 0 0.75 -0.99921 7 0 1.00 -1.52307 8 1 0.00 -0.49070 8 1 0.25 -0.53663 8 1 0.50 -0.65049 9 0 0.00 1.22835 9 0 0.25 1.17369 10 1 0.00 0.83230 11 0 0.00 -0.18064 11 0 0.25 0.08472 11 0 0.50 -1.02813 11 0 0.75 0.33564 12 1 0.00 1.06410 12 1 0.25 1.59505 12 1 0.50 2.05761 12 1 0.75 1.46113 12 1 1.00 1.41302 13 0 0.00 -0.18469 13 0 0.25 -0.18006 13 0 0.50 -0.65446 14 1 0.00 0.05509 15 0 0.00 -1.25072 15 0 0.25 -2.04995 15 0 0.50 -1.36529 15 0 0.75 -2.73013 15 0 1.00 -2.34033 16 1 0.25 0.35311 16 1 0.50 0.51133 16 1 0.75 0.53015 17 0 0.00 1.67120 17 0 0.25 0.61092 18 1 0.00 1.47494 19 0 0.00 -1.00880 19 0 0.25 -0.43964 19 0 0.50 -1.08627 20 1 0.00 1.24841 20 1 0.25 1.04957 20 1 0.50 1.48723 20 1 0.75 2.23080 ; run ; 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; load ; alpha=0.025 ; beta=0.1 ; /* One-sided Type I and Power */ delta=1.178 ; 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,] ; /* Get variance of contrast from fit */ use Soldum1 ; read all into dum1 ; varct=dum1[5,2]**2 ; tau2=nsub*varct ; jackmult=1 ; run jack ; /* calculate jackmult, the inflation factor */ /* 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 ;