About dandar

dandar · ‎08-18-2023

Thanks again

dandar · ‎08-10-2023

@clayt85 Amazing thank you again! I have now familiarised myself with the SAS notation (and this bug). Do you know which is the earliest SAS version that corrects for this?

dandar · ‎08-09-2023

@clayt85 please do not apologise for such a long answer - this is exactly the detail I needed. Thank you so much for taking the time to explain things I greatly appreciate this. The issue of number of subintervals in proc power needing to be high - this makes sense if an approximation is occurring via a piecewise hazard function or something like this (i.e. more pieces is more resolution/granularity) and I suppose I had not yet accepted that, as you correctly point out), sometimes methods just do differ. However 5-6 events can be a big deal in very slow recruiting settings and/or low event rate settings. I suppose "being conservative" and using the largest event number is an option. However your 1-stage (i.e. fixed) sample size design from proc seqdesign does get close so the discrepancy (versus Rpact software in R, and also my manual calculations) is using proc power. I see now I did miss the null hypothesis being tested being -log(HR) and not log(HR), so I suppose I was being lazy and had concluded SAS was "flipping things". However the "sample size not being used" warning from proc seqdesign originating through the duplication of information pertaining to the alternate effect (i.e. using "altref" in the proc design line versus "hazard ratio" in the samplesize line) would probably have baffled me for some time further, so you have saved me a lot of time! Incidentally the same reasons lie behind me using an upper test in proc power - i.e. i realised SAS had flipped the test, and presumably the documentation states it uses -log(HR). On the final subject of the null hazard ratio (my "HR_thresh" macro variable). It is reassuring to know you can obtain the 35 events when using HR_thresh=0.9. Take a look at my output below (I use two versions - one with and without the hazard ratio statement) and they both report a hazard ratio of 0.27 even though both code snippets use HR=0.3. The log(HR/HR_thresh) = log(HR)-log(HR_thresh) should equal -1.0986 but is being reported as -1.20397. The log(HR_thresh) is correctly reported as -0.10536. The log(HR) is incorrectly reported as -1.30933 when it should be -1.2039. Maybe my code is wrong? Can you please post your code - i.e. you may have corrected mine? I am using SAS Studio release 3.81 (Enterprise Edition) SAS release 9.04.01M7P08062020 running on a Linux LIN X64 platform. %let Lambda_cnt=0.01005; %let AT = 1; %let FUT=1; %let HR=0.3; %let HR=0.3; %let alpha=0.025; %let beta=0.1; %let HR_thresh = 0.9; title "Gives HR=0.27: altref specified"; proc seqdesign altref=%sysfunc(log(&HR.)) plots=(none) ; fixed: design nstages=1 alpha=0.025 beta=&beta. alt=upper ; samplesize model=twosamplesurvival (nullhazard=%sysevalf(&Lambda_cnt.*&HR_thresh.) &Lambda_cnt. acctime=&AT. accrual=UNIFORM foltime=&FUT.); run; title "Gives HR=0.27: hazard ratio specified"; proc seqdesign plots=(none) ; fixed: design nstages=1 alpha=0.025 beta=&beta. alt=upper ; samplesize model=twosamplesurvival (nullhazard=%sysevalf(&Lambda_cnt.*&HR_thresh.) &Lambda_cnt. hazardratio=&HR. acctime=&AT. accrual=UNIFORM foltime=&FUT.); run; The output (same for both);

dandar · ‎08-07-2023

I am trying to perform a 1-sided group sequential design using proc seqdesign. I am not too familiar with group sequential designs but believe the issues I am experiencing are not knowing how to "drive the software" rather than GSD-related (I may be wrong however!). In advance I appreciate your time looking at this. The settings are an assumed hazard ratio of 0.3 (<1 favours active treatment) against a threshold of 1 (i.e. testing for superiority at the smallest possible HR of 1). I want 90% power and 2.5% 1-sided type 1 errors. The control arm response rate is 1% events observed in 1 year, and for subject numbers I have chosen a followup time of 1 year and an accrual period of 1 year. The subject numbers is not the primary question, but rather the function "twosamplesurvival" in proc seqdesign uses these inputs. So take these as examples only. Using the following proc power code I get event numbers required of 28 (I "manually" set "numbersubintervals=2" to get 29 to agree with R output of 29 events - increasing this increases event numbers, but this discrepancy is not my primary question - however any thoughts would be welcome). Later I am interested in testing against HR<1 but proc power only seems to test against HR=1 (again any thoughts welcome). %let Lambda_cnt=0.01005; %let AT = 1; %let FUT=1; %let HR=0.3; %let alpha=0.025; %let beta=0.1; proc power; twosamplesurvival test=logrank refsurvexphazard = &Lambda_cnt. hazardratio = &HR. accrualtime = %sysevalf(&AT.) followuptime = %sysevalf(&FUT.) groupweights = (1 1) eventstotal=. alpha=&alpha power = %sysevalf(1-&beta.) sides=U nsubinterval=2; run; The following two code snippets are exactly the same except the first uses "alt=lower" since I want to do a lower 1-sided test. however I get a warning about the sample size procedure not being used. The second snippet runs fine but with an upper boundary with an assumed effect size switched to being positive (the drift parameter in the first is negative as expected, but the second is simply switched to being positive). So my first question is why can I not force SAS to give me a lower test? It appears the results are fine (just "flipped" to be an upper test), but I find this annoying and also unsettling (i.e I am not at all sure I am using the right code). I have been looking at the help documents for a long time and if I have missed a simple option I would be amazed (although I have been known to miss these things!). %let HR_thresh = 1; title "alt=lower: Sample size statement not working"; proc seqdesign altref=%sysfunc(log(&HR.)) stopprob errspend plots=boundary plots=power plots=ASN plots=errspend; ods output Boundary=BI SampleSizeSummary=SSS_D stopprob=stoprob samplesizesummary=sss errspend=errspend; OBrienFleming: design nstages=2 method=obf alpha=&alpha. beta=&beta. alt=lower stop=reject ; samplesize model=twosamplesurvival (nullhazard=%sysevalf(&Lambda_cnt.*&HR_thresh.) &Lambda_cnt. hazardratio=%sysevalf(&HR.) acctime=&AT. accrual=UNIFORM foltime=&FUT.); run; title "alt=upper: Sample size statement working"; proc seqdesign altref=%sysfunc(log(&HR.)) stopprob errspend plots=boundary plots=power plots=ASN plots=errspend; ods output Boundary=BI SampleSizeSummary=SSS_D stopprob=stoprob samplesizesummary=sss errspend=errspend; OBrienFleming: design nstages=2 method=obf alpha=0.025 beta=&beta. alt=upper stop=reject ; samplesize model=twosamplesurvival (nullhazard=%sysevalf(&Lambda_cnt.*&HR_thresh.) &Lambda_cnt. hazardratio=%sysevalf(&HR.) acctime=&AT. accrual=UNIFORM foltime=&FUT.); run; The following is an image of the error for the first code snippet; My second question relates to the threshold. If I set HR_threshold to 0.9 my R sample sizing tells me I need 35 events now rather than 29 due to the more stringent test. But SAS output event numbers stay at 29 (even though a new line appears accounting for this change - circled in red). So what am I doing wrong here?

dandar · ‎01-09-2022

Thank you so much, I will give this a try. I thought I was going mad: interactive mode runs fine (with multiple by-processing errors) but batch mode fails. I'll report back.

dandar · ‎01-09-2022

I had this problem and will try the solutions stated here. However I often use by group processing fitting models on each level (bootstrapping). These often fail in low sample sizes/events and I do wish to simply "run on". Interactive SAS does this but batch does not by default it seems. So I am comfortable with these errors (I capture convergence diagnostics from the sas procedure and filter out non-converging models later). It is best to do a try/catch I suppose but I do not know how, and anyway seems pointless since the base procedure is very happy to just "run on".

dandar · ‎12-10-2020

That indeed did the trick! As I stated earlier, I always thought lsmeans did not do what I needed. Thanks! I am interested in finding out if the NLmixed method is equivalent to what I was going to try in proc IML, but for now I am happy this method works. The code was easy to implement (but time-consuming since I have 35 fixed effects). I chose to supply NLestimate macro with the fixed effects estimates and fixed effects covariance matrix estimate separately that I saved from proc glimmix (as opposed to using the "store" function). I then created a "score" dataset which seems analogous to the "OM" dataset in lsmeans which serves the purpose of supplying all the covariate values need that enter the function expression (I will have multiple rows for multiple comparisons of active doses to placebo). Thus I ended up with; exp(b_c1*x_c1 +...b_cn_p*x_cn_p + b_trt_1*x_trt_1 + ... + b_trt_q*x_trt_q) - exp(b_c1*x_c1 +...b_cn_p*x_cn_p + b_trt_1*x_trt2_1 + ... + b_trt_q*x2_trt2_q) = exp(b_c1*x_c1 +...b_cn_p*x_cn_p) * { exp(b_trt_1*x_trt_1 + ... + b_trt_q*x_trt_q ) - exp(b_trt2_1*x_trt2_1 + ... + b_trt2_q*x_trt2_q ) } where "c" denotes effects common to both terms, and x_trt_1=0 x_trt_2 = 1 ...... x_trt_q = 0 (defines first treatment group of choice) x_trt2_1=1 x_trt2_2 = 0 ...... x_trt2_q = -1 (defines second treatment group of choice) so that the code above defines the mean for treatment group 2 - mean for treatment group q, at the values defined in the common covariates. By creating more rows in the score dataset but changing just the above treatment group dummy variables, each row generates an estimate of an active group - placebo for the same covariates. /*proc glimixx code with manually coded estimate statements The NLmixed macro is subsequently used to estimate differences between the active dose means and the placebo means */ ods output estimates=glmm_p_est covb=glmm_p_FE_cov parameterestimates=glmm_p_FE_est ; proc glimmix data=Model_dat_poiss order=formatted plots=all method=laplace; class usubjid trtp avisitn randstra; model aval = ady trtp avisitn trtp*avisitn base base*ady randstra / dist=poisson link=log solution cl covb ; random int exposure/ subject=usubjid type=UN g gcorr ; output outpred=glmm_p_out pred(blup ilink)=predp lcl(blup ilink)=lowerp ucl(blup ilink)=upperp pred(noblup ilink)=predm lcl(noblup ilink)=lowerm ucl(noblup ilink)=upperm ; covtest glm/ cl; lsmeans trtp / e ilink cl ; estimate 'Severe base: PBO week 1' intercept 1 randstra 1 0 base 21.75 ady 7 base*ady 152.25 trtp 0 0 0 0 1 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0/ilink cl ; estimate 'Severe base: PBO week 4' intercept 1 randstra 1 0 base 21.75 ady 28 base*ady 609 trtp 0 0 0 0 1 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0/ilink cl ; estimate 'Severe base: PBO week 8' intercept 1 randstra 1 0 base 21.75 ady 56 base*ady 1218 trtp 0 0 0 0 1 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0/ilink cl ; estimate 'Severe base: PBO week 12' intercept 1 randstra 1 0 base 21.75 ady 84 base*ady 1827 trtp 0 0 0 0 1 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1/ilink cl ; estimate 'Severe base: 37.5 week 1' intercept 1 randstra 1 0 base 21.75 ady 7 base*ady 152.25 trtp 0 0 1 0 0 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Severe base: 37.5 week 4' intercept 1 randstra 1 0 base 21.75 ady 28 base*ady 609 trtp 0 0 1 0 0 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Severe base: 37.5 week 8' intercept 1 randstra 1 0 base 21.75 ady 56 base*ady 1218 trtp 0 0 1 0 0 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Severe base: 37.5 week 12' intercept 1 randstra 1 0 base 21.75 ady 84 base*ady 1827 trtp 0 0 1 0 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Severe base: 75 week 1' intercept 1 randstra 1 0 base 21.75 ady 7 base*ady 152.25 trtp 0 0 0 1 0 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0/ ilink cl ; estimate 'Severe base: 75 week 4' intercept 1 randstra 1 0 base 21.75 ady 28 base*ady 609 trtp 0 0 0 1 0 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0/ ilink cl ; estimate 'Severe base: 75 week 8' intercept 1 randstra 1 0 base 21.75 ady 56 base*ady 1218 trtp 0 0 0 1 0 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0/ ilink cl ; estimate 'Severe base: 75 week 12' intercept 1 randstra 1 0 base 21.75 ady 84 base*ady 1827 trtp 0 0 0 1 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0/ ilink cl ; estimate 'Severe base: 150 week 1' intercept 1 randstra 1 0 base 21.75 ady 7 base*ady 152.25 trtp 1 0 0 0 0 avisitn 1 0 0 0 trtp*avisitn 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Severe base: 150 week 4' intercept 1 randstra 1 0 base 21.75 ady 28 base*ady 609 trtp 1 0 0 0 0 avisitn 0 1 0 0 trtp*avisitn 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Severe base: 150 week 8' intercept 1 randstra 1 0 base 21.75 ady 56 base*ady 1218 trtp 1 0 0 0 0 avisitn 0 0 1 0 trtp*avisitn 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Severe base: 150 week 12' intercept 1 randstra 1 0 base 21.75 ady 84 base*ady 1827 trtp 1 0 0 0 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Severe base: 300 week 1' intercept 1 randstra 1 0 base 21.75 ady 7 base*ady 152.25 trtp 0 1 0 0 0 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Severe base: 300 week 4' intercept 1 randstra 1 0 base 21.75 ady 28 base*ady 609 trtp 0 1 0 0 0 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Severe base: 300 week 8' intercept 1 randstra 1 0 base 21.75 ady 56 base*ady 1218 trtp 0 1 0 0 0 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Severe base: 300 week 12' intercept 1 randstra 1 0 base 21.75 ady 84 base*ady 1827 trtp 0 1 0 0 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: PBO week 1' intercept 1 randstra 1 0 base 15 ady 7 base*ady 105 trtp 0 0 0 0 1 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0/ilink cl ; estimate 'Moderate base: PBO week 4' intercept 1 randstra 1 0 base 15 ady 28 base*ady 420 trtp 0 0 0 0 1 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0/ilink cl ; estimate 'Moderate base: PBO week 8' intercept 1 randstra 1 0 base 15 ady 56 base*ady 840 trtp 0 0 0 0 1 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0/ilink cl ; estimate 'Moderate base: PBO week 12' intercept 1 randstra 1 0 base 15 ady 84 base*ady 1260 trtp 0 0 0 0 1 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1/ilink cl ; estimate 'Moderate base: 37.5 week 1' intercept 1 randstra 1 0 base 15 ady 7 base*ady 105 trtp 0 0 1 0 0 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: 37.5 week 4' intercept 1 randstra 1 0 base 15 ady 28 base*ady 420 trtp 0 0 1 0 0 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: 37.5 week 8' intercept 1 randstra 1 0 base 15 ady 56 base*ady 840 trtp 0 0 1 0 0 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: 37.5 week 12' intercept 1 randstra 1 0 base 15 ady 84 base*ady 1260 trtp 0 0 1 0 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: 75 week 1' intercept 1 randstra 1 0 base 15 ady 7 base*ady 105 trtp 0 0 0 1 0 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0/ ilink cl ; estimate 'Moderate base: 75 week 4' intercept 1 randstra 1 0 base 15 ady 28 base*ady 420 trtp 0 0 0 1 0 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0/ ilink cl ; estimate 'Moderate base: 75 week 8' intercept 1 randstra 1 0 base 15 ady 56 base*ady 840 trtp 0 0 0 1 0 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0/ ilink cl ; estimate 'Moderate base: 75 week 12' intercept 1 randstra 1 0 base 15 ady 84 base*ady 1260 trtp 0 0 0 1 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0/ilink cl ; estimate 'Moderate base: 150 week 1' intercept 1 randstra 1 0 base 15 ady 7 base*ady 105 trtp 1 0 0 0 0 avisitn 1 0 0 0 trtp*avisitn 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: 150 week 4' intercept 1 randstra 1 0 base 15 ady 28 base*ady 420 trtp 1 0 0 0 0 avisitn 0 1 0 0 trtp*avisitn 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: 150 week 8' intercept 1 randstra 1 0 base 15 ady 56 base*ady 840 trtp 1 0 0 0 0 avisitn 0 0 1 0 trtp*avisitn 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: 150 week 12' intercept 1 randstra 1 0 base 15 ady 84 base*ady 1260 trtp 1 0 0 0 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: 300 week 1' intercept 1 randstra 1 0 base 15 ady 7 base*ady 105 trtp 0 1 0 0 0 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: 300 week 4' intercept 1 randstra 1 0 base 15 ady 28 base*ady 420 trtp 0 1 0 0 0 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: 300 week 8' intercept 1 randstra 1 0 base 15 ady 56 base*ady 840 trtp 0 1 0 0 0 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Moderate base: 300 week 12' intercept 1 randstra 1 0 base 15 ady 84 base*ady 1260 trtp 0 1 0 0 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: PBO week 1' intercept 1 randstra 1 0 base 6.25 ady 7 base*ady 43.75 trtp 0 0 0 0 1 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0/ilink cl ; estimate 'Low base: PBO week 4' intercept 1 randstra 1 0 base 6.25 ady 28 base*ady 175 trtp 0 0 0 0 1 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0/ilink cl ; estimate 'Low base: PBO week 8' intercept 1 randstra 1 0 base 6.25 ady 56 base*ady 350 trtp 0 0 0 0 1 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0/ilink cl ; estimate 'Low base: PBO week 12' intercept 1 randstra 1 0 base 6.25 ady 84 base*ady 525 trtp 0 0 0 0 1 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1/ilink cl ; estimate 'Low base: 37.5 week 1' intercept 1 randstra 1 0 base 6.25 ady 7 base*ady 43.76 trtp 0 0 1 0 0 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 37.5 week 4' intercept 1 randstra 1 0 base 6.25 ady 28 base*ady 175 trtp 0 0 1 0 0 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 37.5 week 8' intercept 1 randstra 1 0 base 6.25 ady 56 base*ady 350 trtp 0 0 1 0 0 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 37.5 week 12' intercept 1 randstra 1 0 base 6.25 ady 84 base*ady 525 trtp 0 0 1 0 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 75 week 1' intercept 1 randstra 1 0 base 6.25 ady 7 base*ady 43.76 trtp 0 0 0 1 0 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 75 week 4' intercept 1 randstra 1 0 base 6.25 ady 28 base*ady 175 trtp 0 0 0 1 0 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 75 week 8' intercept 1 randstra 1 0 base 6.25 ady 56 base*ady 350 trtp 0 0 0 1 0 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0/ilink cl ; estimate 'Low base: 75 week 12' intercept 1 randstra 1 0 base 6.25 ady 84 base*ady 525 trtp 0 0 0 1 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0/ilink cl ; estimate 'Low base: 150 week 1' intercept 1 randstra 1 0 base 6.25 ady 7 base*ady 43.76 trtp 1 0 0 0 0 avisitn 1 0 0 0 trtp*avisitn 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 150 week 4' intercept 1 randstra 1 0 base 6.25 ady 28 base*ady 175 trtp 1 0 0 0 0 avisitn 0 1 0 0 trtp*avisitn 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 150 week 8' intercept 1 randstra 1 0 base 6.25 ady 56 base*ady 350 trtp 1 0 0 0 0 avisitn 0 0 1 0 trtp*avisitn 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 150 week 12' intercept 1 randstra 1 0 base 6.25 ady 84 base*ady 525 trtp 1 0 0 0 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 300 week 1' intercept 1 randstra 1 0 base 6.25 ady 7 base*ady 43.76 trtp 0 1 0 0 0 avisitn 1 0 0 0 trtp*avisitn 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 300 week 4' intercept 1 randstra 1 0 base 6.25 ady 28 base*ady 175 trtp 0 1 0 0 0 avisitn 0 1 0 0 trtp*avisitn 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 300 week 8' intercept 1 randstra 1 0 base 6.25 ady 56 base*ady 350 trtp 0 1 0 0 0 avisitn 0 0 1 0 trtp*avisitn 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; estimate 'Low base: 300 week 12' intercept 1 randstra 1 0 base 6.25 ady 84 base*ady 525 trtp 0 1 0 0 0 avisitn 0 0 0 1 trtp*avisitn 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0/ilink cl ; run; /* create just one row to generate a difference of treatment group dose 37.5 - placebo*/ data est_cov_dat; length base_type $ 15. ady 4 dose 4 intercept 4 ady 4 trt_150 4 trt_300 4 trt_375 4 trt_75 4 trt_pbo 4 avisitn_3 4 avisitn_6 4 avisitn_10 4 avisitn_14 4 trt_150_avisitn_3 4 trt_150_avisitn_6 4 trt_150_avisitn_10 4 trt_150_avisitn_14 4 trt_300_avisitn_3 4 trt_300_avisitn_6 4 trt_300_avisitn_10 4 trt_300_avisitn_14 4 trt_375_avisitn_3 4 trt_375_avisitn_6 4 trt_375_avisitn_10 4 trt_375_avisitn_14 4 trt_75_avisitn_3 4 trt_75_avisitn_6 4 trt_75_avisitn_10 4 trt_75_avisitn_14 4 trt_pbo_avisitn_3 4 trt_pbo_avisitn_6 4 trt_pbo_avisitn_10 4 trt_pbo_avisitn_14 4 trt2_150 4 trt_300 4 trt2_375 4 trt_75 4 trt2_pbo 4 trt2_150_avisitn_3 4 trt2_150_avisitn_6 4 trt2_150_avisitn_10 4 trt2_150_avisitn_14 4 trt2_300_avisitn_3 4 trt2_300_avisitn_6 4 trt2_300_avisitn_10 4 trt2_300_avisitn_14 4 trt2_375_avisitn_3 4 trt2_375_avisitn_6 4 trt2_375_avisitn_10 4 trt2_375_avisitn_14 4 trt2_75_avisitn_3 4 trt2_75_avisitn_6 4 trt2_75_avisitn_10 4 trt2_75_avisitn_14 4 trt2_pbo_avisitn_3 4 trt2_pbo_avisitn_6 4 trt2_pbo_avisitn_10 4 trt2_pbo_avisitn_14 4 base 4 ady_base 4 randstra_comp 4 randstra_noncomp 4 ; base_type = "Severe"; ady = 7; dose = 37.5; intercept = 1; ady = 7; trt_150 = 0; trt_300 = 0; trt_375 = 1; trt_75 = 0; trt_pbo = 0; trt2_150 = 0; trt2_300 = 0; trt2_375 = 0; trt2_75 = 0; trt2_pbo = -1; avisitn_3 = 1; avisitn_6 = 0; avisitn_10 = 0; avisitn_14 = 0; trt_150_avisitn_3 = 0; trt_150_avisitn_6 = 0; trt_150_avisitn_10 = 0; trt_150_avisitn_14 = 0; trt_300_avisitn_3 = 0; trt_300_avisitn_6 = 0; trt_300_avisitn_10 = 0; trt_300_avisitn_14 = 0; trt_375_avisitn_3 = 1; trt_375_avisitn_6 = 0; trt_375_avisitn_10 = 0; trt_375_avisitn_14 = 0; trt_75_avisitn_3 = 0; trt_75_avisitn_6 = 0; trt_75_avisitn_10 = 0; trt_75_avisitn_14 = 0; trt_pbo_avisitn_3 = 0; trt_pbo_avisitn_6 = 0; trt_pbo_avisitn_10 = 0; trt_pbo_avisitn_14 = 0; trt2_150_avisitn_3 = 0; trt2_150_avisitn_6 = 0; trt2_150_avisitn_10 = 0; trt2_150_avisitn_14 = 0; trt2_300_avisitn_3 = 0; trt2_300_avisitn_6 = 0; trt2_300_avisitn_10 = 0; trt2_300_avisitn_14 = 0; trt2_375_avisitn_3 = 0; trt2_375_avisitn_6 = 0; trt2_375_avisitn_10 = 0; trt2_375_avisitn_14 = 0; trt2_75_avisitn_3 = 0; trt2_75_avisitn_6 = 0; trt2_75_avisitn_10 = 0; trt2_75_avisitn_14 = 0; trt2_pbo_avisitn_3 = -1; trt2_pbo_avisitn_6 = 0; trt2_pbo_avisitn_10 = 0; trt2_pbo_avisitn_14 = 0; base = 21.75; ady_base = 152.25; randstra_comp=1; randstra_noncomp=0; output; run; /*run the NLestimate macro*/ %NLEstimate(inest=glmm_p_FE_est,incovb=glmm_p_FE_cov_num, label=Rate Difference, f=exp(b_p1*intercept + b_p2*ady + b_p8*avisitn_3 + b_p9*avisitn_6 + b_p10*avisitn_10 + b_p11*avisitn_14 + b_p32*base + b_p33*ady_base + b_p34*randstra_comp + b_p35*randstra_noncomp) * (exp( b_p3*trt_150 + b_p4*trt_300 + b_p5*trt_375 + b_p6*trt_75 + b_p7*trt_pbo + b_p12*trt_150_avisitn_3 + b_p13*trt_150_avisitn_6 + b_p14*trt_150_avisitn_10 + b_p15*trt_150_avisitn_14 + b_p16*trt_300_avisitn_3 + b_p17*trt_300_avisitn_6 + b_p18*trt_300_avisitn_10 + b_p19*trt_300_avisitn_14 + b_p20*trt_375_avisitn_3 + b_p21*trt_375_avisitn_6 + b_p22*trt_375_avisitn_10 + b_p23*trt_375_avisitn_14 + b_p24*trt_75_avisitn_3 + b_p25*trt_75_avisitn_6 + b_p26*trt_75_avisitn_10 + b_p27*trt_75_avisitn_14 + b_p28*trt_pbo_avisitn_3 + b_p29*trt_pbo_avisitn_6 + b_p30*trt_pbo_avisitn_10 + b_p31*trt_pbo_avisitn_14 )- exp(b_p3*trt2_150 + b_p4*trt2_300 + b_p5*trt2_375 + b_p6*trt2_75 + b_p7*trt2_pbo + b_p12*trt2_150_avisitn_3 + b_p13*trt2_150_avisitn_6 + b_p14*trt2_150_avisitn_10 + b_p15*trt2_150_avisitn_14 + b_p16*trt2_300_avisitn_3 + b_p17*trt2_300_avisitn_6 + b_p18*trt2_300_avisitn_10 + b_p19*trt2_300_avisitn_14 + b_p20*trt2_375_avisitn_3 + b_p21*trt2_375_avisitn_6 + b_p22*trt2_375_avisitn_10 + b_p23*trt2_375_avisitn_14 + b_p24*trt2_75_avisitn_3 + b_p25*trt2_75_avisitn_6 + b_p26*trt2_75_avisitn_10 + b_p27*trt2_75_avisitn_14 + b_p28*trt2_pbo_avisitn_3 + b_p29*trt2_pbo_avisitn_6 + b_p30*trt2_pbo_avisitn_10 + b_p31*trt2_pbo_avisitn_14)), score=est_cov_dat,title=Difference means,outscore=score_out); proc print data=score_out; var base_type ady dose pred StdErrPred lower upper ChiSq pr; run; Output for dose 37.5 - placebo for severe baseline subjects at day 7; Group estimates from glimmix - we can see the correct difference in means was calculated

dandar · ‎12-10-2020

Thanks so much for your help.

dandar · ‎12-10-2020

Ahh! So lsmeans with the "diff" option computes the difference on the natural scale with the correct variance - i.e. it calculates exp(eta1)-exp(eta2) rather than exp[eta1-eta2]? I think I have had misconceptions about lsmeans then for a long time! I always tend to use the "estimate" statement, and this calculates exp[eta1-eta2] rather like lsmestimate does. Thus I convinced myself lsmeans would do the same (or I forgot it can do better).

dandar · ‎12-10-2020

Thanks for your quick response. Sorry yes my data are Poisson distributed. I know the common quantity computed are ratios of group means (so the other covariate effects cancel out under the log link), but clinical folks like to see differences. For this reason I have for a long time wanted to be able to at least have the ability to compute these differences and the variance so I can easily construct an asymptotic Normal CI on it. You think this can be done with a custom variance function?

dandar · ‎12-10-2020

I am interested in being able to calculate differences in means of groups (i.e. treatment - control) on the natural scale from within Proc Glimmix. For instance the snippet below from the help document for the "lsmestimate" statement suggests only the difference on the linear predictor scale can be computed. I am in fact fitting a repeated measures model to count data so am interested in the log link, however the principle is the same. I could work out the variance of the difference on the natural scale, and this would (I think) be a function of the covariance matrix for the fixed effects. I could I suppose therefore output this estimated covariance matrix, and fixed effects estimates, read it all into IML, and compute the variance of this difference "manually". However....if there is an easier way from within Proc Glimmix? Any help would be greatly appreciated.

Online Status	Offline
Date Last Visited	‎08-30-2023 04:38 PM

Re: Problems performing 1-sided group sequential design in proc seqdes...

Re: Problems performing 1-sided group sequential design in proc seqdes...

Re: Problems performing 1-sided group sequential design in proc seqdes...

Problems performing 1-sided group sequential design in proc seqdesign

Re: Problem with sas batch mode

Re: Problem with sas batch mode

Re: Calculating variance of differences in Proc Glimmix on the natural...

Re: Calculating variance of differences in Proc Glimmix on the natural...

Re: Calculating variance of differences in Proc Glimmix on the natural...

Re: Calculating variance of differences in Proc Glimmix on the natural...

Re: Problem with sas batch mode

Re: Calculating variance of differences in Proc Glimmix on the natural...

Re: Problems performing 1-sided group sequential design in proc seqdes...

Re: Problems performing 1-sided group sequential design in proc seqdes...

Re: Problems performing 1-sided group sequential design in proc seqdes...

Problems performing 1-sided group sequential design in proc seqdesign

Re: Problem with sas batch mode

Re: Problem with sas batch mode

Re: Calculating variance of differences in Proc Glimmix on the natural...

Re: Calculating variance of differences in Proc Glimmix on the natural...

Re: Calculating variance of differences in Proc Glimmix on the natural...

Re: Calculating variance of differences in Proc Glimmix on the natural...

Calculating variance of differences in Proc Glimmix on the natural sca...