95% CI of relative change

[Lindsy]

Hi,
I have a problem related to the calculation of 95% CI of relative change.
I ran an interrupted time series using segmented linear regression model and plan to show the absolute and relative change at one point of time however, I don’t know how to run 95% CI of relative change.
I found one SAS code. They used bootstrap and delta method to calculate this 95% CI of relative change. But I’m not good at it. Please suggest what I should do?
Thank you in advance!

[Dimitriy V. Masterov]

Can you use nlcom to calculate the relative change? Like nlcom rel_delta: (_b[after]-_b[before])/(_b[before])?

It would help if you could post code with a fake dataset.

[Lindsy]

Thank you for your answer. Yes, I can calculate relative change using nlcom.
However, how to calculate 95%CI around that result.
For, SAS code. Two methods were used to calculate bootstrap and delta methods use for 95% CI calculation

options ls=80; libname dat 'C:\Work Projects\Timeseries_CI\Final texts numbers pgms and outputs'; /*---------------------------; -----Total Prescription------; --------------------------- ;*/ data drug1;set dat.red_war2; month=_n_; intervention=(month>=36); monthafintervention=max(0,month-35); run; data drug2;set drug1; if 39>=month >=36 then delete; if month >= 36 then month=month-4; monthafintervention=max(month-35,0); run; proc autoreg data=drug2; model war_other5=month intervention monthafintervention/backstep slstay=0.05 nlag=12 method=ml covb; model war_other5=intervention monthafintervention/backstep slstay=0.05 nlag=12 method=ml covb; run; *---------------; *Delta method; *---------------; *Absolute and relative change and 95% CI for level change; %diffrr(drug2,drg_lv,war_other5,intervention monthafintervention,12,1 0 0,1 1 0,0.05); *Absolute and relative change with 95% CI six months after the alert; %diffrr(drug2,drg_sixm,war_other5,intervention monthafintervention,12,1 0 0,1 1 6,0.05); *---------------------; *Bootstrapping method; *---------------------; options nonotes nosyntaxcheck nosymbolgen; %boots(3263,0,-311.3794,-21.2870,52, 25993,%str(1),10000,36,12, time intervention1 timeafterintervention1); data meanvar; set simulationresults(firstobs=2); lvrel=intervention1/intercept; *slrel=timeafterintervention1/time; sixmabs=intervention1+6*timeafterintervention1; sixmrel=(intervention1+6*timeafterintervention1)/(intercept+41*time); run; %macro nm(vr); proc univariate data=meanvar alpha=0.05; var &vr; output out=&vr pctlpre=col pctlpts=2.5,50,97.5 mean=me ; run; proc print;run; %mend nm; *Use median is better; %nm(lvrel);%nm(sixmrel); /*-----------------------------------; -----Acetaminophen Prescription------; ------------------------------------;*/ *---------------; *Delta method; *---------------; data drug1;set dat.red_war2; month=_n_; intervention=(month>=36); monthafintervention=max(0,month-35); run; data drug2;set drug1; if 39>=month >=36 then delete; if month >= 36 then month=month-4; monthafintervention=max(month-35,0); run; proc autoreg data=drug2; model war_apap=month intervention monthafintervention/backstep slstay=0.05 nlag=12 method=ml covb; run; *Absolute and relative change and 95% CI for level change; %diffrr(drug2,drg_lv,war_apap,month intervention monthafintervention,12,1 0 0 0,1 0 1 0,0.05); *Absolute and relative change and 95% CI for trend change; %diffrr(drug2,drg_slope,war_apap,month intervention monthafintervention,12,0 1 0 0,0 1 0 1,0.05); *Absolute and relative change with 95% CI six months after the alert; %diffrr(drug2,drg_sixm,war_apap,month intervention monthafintervention,12,1 41 0 0,1 41 1 6,0.05); *---------------------; *Bootstrapping method; *---------------------; options nonotes nosyntaxcheck nosymbolgen; %boots(2321,8.2995,-310.9912,-22.1976,52, 24017,%str(1),10000,36,12, time intervention1 timeafterintervention1); data meanvar; set simulationresults(firstobs=2); lvrel=intervention1/intercept; slrel=timeafterintervention1/time; sixmabs=intervention1+6*timeafterintervention1; sixmrel=(intervention1+6*timeafterintervention1)/(intercept+41*time); run; %macro nm(vr); proc univariate data=meanvar alpha=0.05; var &vr; output out=&vr pctlpre=col pctlpts=2.5,50,97.5 mean=me ; run; proc print;run; %mend nm; *Use median is better; %nm(lvrel);%nm(slrel);%nm(sixmrel); ***SAS code for bootstrap
/* The macro %BOOTS constructs confidence intervals around the absolute and relative effect estimates from segmented linear regressions using bootstrapping method. We simulate a time series based on the final selected model at each iteration. We route the output to a pre-designated file so that we can exclude iterations where the maximum likelihood method does not converge. We rerun the segmented time series regression to obtain estimates for each valid iteration. We summarize the distribution of all estimates from the simulations and construct the confidence interval from data set simulation results with flag delete=0. 1. Identify a temporary file to which the SAS output will be routed by changing the path in the first statement. (For example: filename routed 'c:\zzzz.txt' 2. Specify the following required parameters for the macro: Baseint - the estimated baseline intercept from the final model Baseslope - the estimated baseline slope from the final model Interaint - the estimated level change comparing the intervention period to the baseline from the final model Interaslope - the estimated trend change comparing the intervention period to the baseline from the final model Nperiod - the number of total time points Sigmas - the estimated variance from the final model Ar - the estimated autocorrelation terms from the final model Iter - the number of simulation iterations Intma - the value of the starting time point of the first policy Modar - the maximum allowable autocorrelations will be used in each iteration of segmented regression desired by the investigators Fieldb - the list of significant terms from the final model */ data logall;intercept=0;time=0; intervention1=0;timeafterintervention1=0; delete=1; run; %macro BOOTS(baseint,baseslope,interaint,interaslope, nperiod,sigmas,ar,iter,intma,modar,fieldb); %macro numa;%do i=1 %to &iter;logd&i%end;%mend numa; data xbeta(drop=beta0-beta3); beta0=&baseint;beta1=&baseslope;beta2=&interaint;b eta3=&interaslope; do time=1 to &nperiod; intervention1=(time >=&intma); timeafterintervention1=max(0,time-&intma+1); mu2=beta0+beta1*time+beta2*intervention1+beta3*tim eafterintervention1; output; end; run; %do k=1 %to &iter; proc iml; reset noname; use xbeta var {mu2}; read all into z2; start armasim(y,n,phi,theta,seed,sigma_2); p = ncol(phi)-1;q = ncol(theta)-1; y = normal(j(1,n+q,seed));y = sqrt(sigma_2)*y; /*---- Pure MA or white noise ----*/ if p=0 then y=product(theta,y)[,(q+1)n+q)]; else do; /*---- Pure AR or ARMA ----*/ /*---- Get the autocovariance function ----*/ call armacov(gamma,cov,ma,phi,theta,p); if gamma[1]<0 then do; print 'ARMA parameters not stable.'; print 'Execution terminating.'; stop; end; /*---- Form covariance matrix ----*/ gamma=toeplitz(gamma); /*---- Generate covariance between initial y and ----*/ /*---- initial innovations ----*/ if q>0 then do; psi = ratio(phi,theta,q); psi = hankel( psi[,-( (-q)-1) ) ] ); m = max(1,(q-p+1)); psi = psi[ -((-q)-m)) ,]; if p>q then psi = j(p-q,q,0)//psi; gamma = ( gamma||psi ) // ( psi`||i(q) ); end; /*---- Use Cholesky root to get startup values ----*/ gamma = root(gamma); startup = y[,1p+q)]*gamma; e = y[,(p+q+1)n+q)]; /*---- Generate MA part ----*/ if q>0 then do; e = startup[,(p+1)p+q)]||e; e = product(theta,e)[,(q+1)n+q-p)]; end; y = startup[,1:p]; phi1 = phi[, -(-(p+1)-2)) ]`; /*---- Use difference equation to generate ----*/ /*---- remaining values ----*/ do ii=1 to n-p; y = y||( e[,ii]-y[,iiii+p-1)]*phi1 ); end; end; y = y`; finish; /*---- ARMASIM ----*/ seed=%eval(&k*&iter); sample=%eval(&nperiod+200); sam=sample+1-%eval(&nperiod); run armasim(y,sample,{&ar},{1},seed,&sigmas); ya=y[sam:sample]; yf=ya+z2; create ts from yf; append from yf; quit; data dd; merge xbeta ts; /*if &phaseup ne . then do;if &phaseup>=time>=&phaselow then col1=.;end; else do;end;*/ run; filename routed 'c:\zzzz.txt'; proc printto print=routed new; run; proc autoreg data=dd outest=logd(keep=intercept &fieldb); model col1=&fieldb/backstep method=ml nlag=&modar; run; proc printto;run; data test; infile routed pad end=eof; input word $char256.; retain ml mle 0; if index(word,'Maximum Likelihood Estimates') > 0 then ml=1; if index(word,'ERROR: The estimation did not converge') >0 then mle=1; if eof then call symput('ml',compress(ml)) ; if eof then call symput('mle',compress(mle)) ; run; %if &ml=1 and &mle=1 %then %do; data logd;intercept=0;time=0; intervention1=0;timeafterintervention1=0; delete=1; run; proc append base=logall data=logd;run; %end; %else %if &ml=1 and &mle ne 1 %then %do; data logd;set logd;delete=0; proc append base=logall data=logd;run; %end; %else %do; data logd;set logd;delete=0; proc append base=logall data=logd;run; %end; %end; data simulationresults;set logall;run; %mend BOOTS; **SAS for Delta Method/*
The macro %DIFFRR constructs confidence intervals around the absolute and relative effect estimates from segmented linear regressions using Delta method. We run the selected final segmented autoregressive error model, routing the output to a pre-designated file and writing to data sets the parameter estimates and covariance matrix obtained through maximum likelihood estimation, assuming as given the autocorrelations estimated by the procedure. Then we choose the maximum likelihood estimates, if available. Otherwise, we use OLS estimates if maximum likelihood estimation does not converge. One can use Yule-Walker or other types of estimates by changing the corresponding codes. Finally, we use the IML procedure to calculate the expectation and variance of the absolute difference as well as the relative change using the multivariate delta method for the latter. Normal distribution (based on large sample theory) is assumed to construct the final confidence intervals. 1. Identify a temporary file to which the SAS output will be routed by changing the path in the first statement. (For example: filename routed 'c:\zzzz.txt' 2. Specify the following seven required parameters for the macro: In - the name of the time series data set to be used to run the final model Out - the output data set with the information created by the macro Outcome - the dependent variable of the time series regression model Predictors - the independent variables in the final time series regression model Lag - the lag term for the final time series regression model Variablewithout - the values of the variables in the final model at the time point of interest, without the intervention Variablewith - the values of the variables in the final model at the time point of interest, with the intervention */ %macro DIFFRR(in,out,outcome,predictors,lag,coefwithout,c oefwith,siglevel); %macro sort(dsn,ods,by,where,keep); proc sort data=&dsn(keep=&keep) out=&ods; by &by; where &where; run; %mend sort; filename routed 'c:\zzzz.txt'; proc printto print=routed new; run; /*we output the parameter estimates and covariate matrix of the final ML model that assumes the autregressive parameters given*/ proc autoreg data=&in outest=est covout; ods output Autoreg.Model1.FinalModel.CovB=work.afinalcov&in; ods output Autoreg.Model1.FinalModel.ParameterEstimatesGivenA R=work.finalparam&in(keep=estimate variable); ods output Autoreg.Model1.OLSEst.CovB=work.olscov&in(drop=mod el); ods output Autoreg.Model1.OLSEst.ParameterEstimates=work.olsp aram&in(keep=estimate variable); model &outcome=&predictors/ method=ml nlag=&lag backstep dwprob covb; run; proc printto; run; data test; infile routed pad end=eof; input word $char256.; retain ml mle 0; if index(word,'Maximum Likelihood Estimates') > 0 then ml=1; if index(word,'ERROR: The estimation did not converge') >0 then mle=1; if eof then call symput('ml',compress(ml)) ; if eof then call symput('mle',compress(mle)) ; run; %if &ml=1 and &mle ne 1 %then %do; data afinalcov&in;set afinalcov&in(keep=rowname intercept &predictors); j=_n_; run; %sort(afinalcov&in,tfinalcov&in,rowname j); data finalcov&in(where=(substr(rowname,1,2) ne 'AR'));set tfinalcov&in;by rowname;if last.rowname;run; %sort(finalcov&in,finalcov&in,j); data d1;set finalcov&in(drop=j);run; data d2;set finalparam&in;run; %end; %else %do;data d1;set olscov&in;run; data d2;set olsparam&in;run; %end; /*we use the output of the parameter estimates of the final model*/ proc transpose data=d2 out=param(drop=_name_); var estimate; id variable; run; proc iml; use d1; read all into x; use param; read all into est; c={&coefwithout}; d={&coefwith}; dis=d-c; ywithout=c*est`; ywith=d*est`; varwithout=c*x*c`; varwith_dis=dis*x*dis`; cov=c*x*dis`; absdiff=(ywith-ywithout); relchange=(ywith-ywithout)/ywithout; y=j(1,10); varabsdiff=dis*x*dis`; varrelchange=(1/(ywithout**2))*varwith_dis+(absdiff**2/(ywithout**4))*varwithout- 2*(absdiff/(ywithout**3))*cov; hfcit=1-(&siglevel)/2; crt=quantile('NORMAL',hfcit); lowerabsdiff=absdiff-crt*sqrt(varabsdiff); upperabsdiff=absdiff+crt*sqrt(varabsdiff); lowerrelchange=relchange-crt*sqrt(varrelchange); upperrelchange=relchange+crt*sqrt(varrelchange); y[1,1]=ywithout;y[1,2]=ywith; y[1,3]=absdiff;y[1,4]=varabsdiff;y[1,5]=lowerabsdiff;y[1,6]=upperabsdiff; y[1,7]=relchange;y[1,8]=varrelchange;y[1,9]=lowerrelchange;y[1,10]=upperrelchange; create variance from y; append from y; quit; data &out;set variance;rename col1=y1 col2=y2 col3=absdiff col4=varabsdiff col5=lowerabsdiff col6=upperabsdiff col7=relchange col8=varrelchange col9=lowerrelchange col10=upperrelchange; run; %mend DIFFRR;
Please suggest how should I convert these SAS codes for STATA.

[Dimitriy V. Masterov]

nlcom uses the delta method to give you the CI.

If that's not what you want, it may help if you precisely define what a "CI around that result means" using math or words.

I tried to parse that SAS code, but it it formatted too poorly to be legible.