2
$\begingroup$

I am simulating a logistic regression setting (P(Y=1|x). Estimated proportion of patients with a poor feature at time zero is 40%. I am generating 1000 replications per scenario, sample size ranging from 100 to 200 by 10, for odds ratios ranging from 1.5 to 2.5 by 0.1. I generate x (years) from Uniform~(2,6). Relative power increases appropriately for odds ratios from 1.5 to 2.0. However, relative power decreases for odds ratios greater than 2.0. Does anyone have an explanation for this.

enter image description here

data t; 
b0=-0.405;`enter code here`
do or=1.5 to 2.5 by 0.1;
b1=log(or);
do ss=100 to 200 by 10;
    do rep=1 to 1000;
        do i=1 to ss;
        *years=6*ranuni(34958487);
        years=4*ranuni(34958487)+2;
        *years=1*rannor(34958487)+3;
        lp=b0 + b1*years;
        pi=exp(lp)/(1 + exp(lp));
        y=ranbin(45823765,1,pi);
        output;
        end;
    end;
end;
end;
run;

proc sort data=t; by b1 ss rep;run;
ods listing close;
ods select ParameterEstimates(persist);
ods output ParameterEstimates=Estimates;
proc logistic descending data=t;
by b1 ss rep;
model y=years;
run;
ods select all;
ods listing;

data reject; set estimates;
if variable="years";
rej=(probchisq<0.05);
or=exp(b1);
run;
proc sort data=reject; by or ss;run;
options pageno=1 pagesize=max;
proc freq data=reject;
by or;
tables ss*rej/ out=power;
title "Rejection Probability";
run;

data power; set power;
if rej=1;
rejprob=count/1000;
keep or ss rejprob;
run;
$\endgroup$

1 Answer 1

1
$\begingroup$

This is due to ceiling effects. Look at a plot of the success probabilities in some of your simulated data sets when $\beta_1$ is near the higher end of your range, $\log(2.4)$, vs. the lower end, $\log(1.6)$. Toward the higher end, far more of the predicted probabilities get near upper bound of 1. For example, when $\beta_1 = \log(1.6)$, around $15\%$ of predicted probabilities will be above $.95$; when $\beta_1 = \log(2.4)$, that is more like $50\%$. This dynamic makes it harder to determine the precise effect of the predictor, which is reflected in an increased standard error, and thus a smaller power. To get the expected results, you should fix the marginal probability, $P(Y=1)$, across simulations by adjusting $\beta_0$ appropriately.

$\endgroup$
3
  • $\begingroup$ I fixed b0 such that the probability of the event at time zero was 40% in all models. I am confident your diagnosis is correct. However, I do not know how to fix the marginal probability, P(Y=1), across simulations by adjusting b0. Can you expand on this? Thank you. $\endgroup$
    – Dan Ayers
    Commented Feb 21, 2016 at 6:55
  • $\begingroup$ The only thing that comes to mind is that you could calculate $$ P(Y=1) = E_X (P(Y=1|X)) $$ as a function of $\beta_0, \beta_1$. Given the distribution of $X$, that would equal $$ f(\beta_0, \beta_1) = \frac{1}{4} \int_{2}^{6} \frac{1}{1 + e^{-(\beta_0 + \beta_1 x)}} dx $$ I don't think you can do this analytically but, for a fixed value of $\beta_1$ and $P(Y=1)$, you can find the proper value of $\beta_0$ by numerically solving for a root of the equation $$f(\beta_0, \beta_1) - P(Y=1) $$ $\endgroup$
    – linksys
    Commented Feb 22, 2016 at 4:00
  • $\begingroup$ You could do fairly well by setting $P(Y=1)$ at the midpoint of time to a fixed value, which you can do by adjusting the intercept until it's right. $\endgroup$ Commented May 23 at 21:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.