Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Power calculation for interaction effect

    Hi,

    I am trying to do a power calculation for an interaction effect.
    My outcome is incident CVD and my exposure is social asymmetry category (sacat) in 4 groups. I want to see if country group (grouped according to welfare) interacts with social asymmetry and I'm trying to do my power calculation as the numbers are very small for some of the social asymetry groups.
    This is my code:

    capture program drop sim_interaction_sacat
    program sim_interaction_sacat, rclass
    version 16.1
    syntax, n(integer) alpha(0.05)

    //set sample size

    set obs `n'

    // Group ID assignment: 16 total combinations (4 sacat x 4 welfare)
    gen combo = runiform()
    gen group = .

    // Assign each combo based on observed proportions

    replace group = 1 if combo < .24312135 // Sacat 1 / Welfare 1
    replace group = 2 if combo >= .24312135 & combo < .32511415 // Sacat 2 / Welfare 1
    replace group = 3 if combo >= .32511415 & combo < .39286705 // Sacat 3 / Welfare 1
    replace group = 4 if combo >= .39286705 & combo < 0.4535 // Sacat 4 / Welfare 1

    replace group = 5 if combo >= 0.4535 & combo < .560189 // Sacat 1 / Welfare 2
    replace group = 6 if combo >= .560189 & combo < .59755028 // Sacat 2 / Welfare 2
    replace group = 7 if combo >= .59755028 & combo < .62979854 // Sacat 3 / Welfare 2
    replace group = 8 if combo >= .62979854 & combo < 0.6549 // Sacat 4 / Welfare 2

    replace group = 9 if combo >= 0.6549 & combo < .71310498 // Sacat 1 / Welfare 3
    replace group = 10 if combo >= .71310498 & combo < .74304366 // Sacat 2 / Welfare 3
    replace group = 11 if combo >= .74304366 & combo < .75589944 // Sacat 3 / Welfare 3
    replace group = 12 if combo >= .75589944 & combo < .7575 // Sacat 4 / Welfare 3

    replace group = 13 if combo >= .7575 & combo < .8559065 // Sacat 1 / Welfare 4
    replace group = 14 if combo >= .8559065 & combo < .93270625 // Sacat 2 / Welfare 4
    replace group = 15 if combo >= .93270625 & combo < .96908125 // Sacat 3 / Welfare 4
    replace group = 16 if combo >= .96908125 // Sacat 4 / Welfare 4

    // Map to groups
    gen sacat = ceil(mod(group-1,4)+1)
    gen welfare = ceil((group-1)/4 + 1)

    // Assign observed incident CVD proportions to groups
    gen p = .
    replace p = 0.0503 if group == 1
    replace p = 0.0645 if group == 2
    replace p = 0.0504 if group == 3
    replace p = 0.0658 if group == 4

    replace p = 0.0823 if group == 5
    replace p = 0.0924 if group == 6
    replace p = 0.0387 if group == 7
    replace p = 0.1281 if group == 8

    replace p = 0.0672 if group == 9
    replace p = 0.0823 if group == 10
    replace p = 0.0874 if group == 11
    replace p = 0.0808 if group == 12

    replace p = 0.0702 if group == 13
    replace p = 0.1002 if group == 14
    replace p = 0.0698 if group == 15
    replace p = 0.1252 if group == 16

    gen outcome = runiform() < p

    // Run logistic regression with interaction
    logit outcome i.sacat##i.welfare, nolog
    testparm i.sacat#i.welfare

    return scalar reject = (r(p) < `alpha')
    clear
    end

    simulate reject = r(reject), reps(500) seed(12345): ///
    sim_interaction_sacat, n(24031) alpha(0.05)


    But its coming up with this error message:
    . simulate reject = r(reject), reps(500) seed(12345): ///
    > sim_interaction_sacat, n(24031) alpha(0.05)
    invalid syntax
    an error occurred when simulate executed sim_interaction_sacat
    r(197);

    end of do-file

    r(197);

    and I can't work out why!

    I've tried to change the program to use args instead of syntax (args n alpha.... sim_interaction_sacat 24031 0.05)but that also comes up with the same error message

    Thanks very much for any help
    Last edited by Isabel Westlake; 15 Jul 2025, 05:54. Reason: regression

  • #2
    I only read r(197) and the short program at the top of your post, so I won't comment on content. Syntax-wise, you want
    Code:
    syntax, n(integer) [ alpha(real 0.05) ]

    Comment


    • #3
      Ahh sorry has this worked?
      Code:
      capture program drop sim_interaction_sacat
      program sim_interaction_sacat, rclass
      version 16.1
      syntax, n(integer) alpha(0.05)
      
      //set sample size
      
      set obs `n'
      
      // Group ID assignment: 16 total combinations (4 sacat x 4 welfare)
      gen combo = runiform()
      gen group = .
      
      // Assign each combo based on observed proportions
      
      replace group = 1 if combo < .24312135 // Sacat 1 / Welfare 1
      replace group = 2 if combo >= .24312135 & combo < .32511415 // Sacat 2 / Welfare 1
      replace group = 3 if combo >= .32511415 & combo < .39286705 // Sacat 3 / Welfare 1
      replace group = 4 if combo >= .39286705 & combo < 0.4535 // Sacat 4 / Welfare 1
      
      replace group = 5 if combo >= 0.4535 & combo < .560189 // Sacat 1 / Welfare 2
      replace group = 6 if combo >= .560189 & combo < .59755028 // Sacat 2 / Welfare 2
      replace group = 7 if combo >= .59755028 & combo < .62979854 // Sacat 3 / Welfare 2
      replace group = 8 if combo >= .62979854 & combo < 0.6549 // Sacat 4 / Welfare 2
      
      replace group = 9 if combo >= 0.6549 & combo < .71310498 // Sacat 1 / Welfare 3
      replace group = 10 if combo >= .71310498 & combo < .74304366 // Sacat 2 / Welfare 3
      replace group = 11 if combo >= .74304366 & combo < .75589944 // Sacat 3 / Welfare 3
      replace group = 12 if combo >= .75589944 & combo < .7575 // Sacat 4 / Welfare 3
      
      replace group = 13 if combo >= .7575 & combo < .8559065 // Sacat 1 / Welfare 4
      replace group = 14 if combo >= .8559065 & combo < .93270625 // Sacat 2 / Welfare 4
      replace group = 15 if combo >= .93270625 & combo < .96908125 // Sacat 3 / Welfare 4
      replace group = 16 if combo >= .96908125 // Sacat 4 / Welfare 4
      
      // Map to groups
      gen sacat = ceil(mod(group-1,4)+1)
      gen welfare = floor((group-1)/4 + 1)
      
      // Assign observed incident CVD proportions to groups
      gen p = .
      replace p = 0.0503 if group == 1
      replace p = 0.0645 if group == 2
      replace p = 0.0504 if group == 3
      replace p = 0.0658 if group == 4
      
      replace p = 0.0823 if group == 5
      replace p = 0.0924 if group == 6
      replace p = 0.0387 if group == 7
      replace p = 0.1281 if group == 8
      
      replace p = 0.0672 if group == 9
      replace p = 0.0823 if group == 10
      replace p = 0.0874 if group == 11
      replace p = 0.0808 if group == 12
      
      replace p = 0.0702 if group == 13
      replace p = 0.1002 if group == 14
      replace p = 0.0698 if group == 15
      replace p = 0.1252 if group == 16
      
      gen outcome = runiform() < p
      
      // Run logistic regression with interaction
      logit outcome i.sacat##i.welfare, nolog
      testparm i.sacat#i.welfare
      
      return scalar reject = (r(p) < `alpha')
      clear
      end
      
      simulate reject = r(reject), reps(500) seed(12345): ///
      sim_interaction_sacat, n(24031) alpha(0.05)
      Last edited by Isabel Westlake; 15 Jul 2025, 07:09.

      Comment


      • #4
        Originally posted by daniel klein View Post
        I only read r(197) and the short program at the top of your post, so I won't comment on content. Syntax-wise, you want
        Code:
        syntax, n(integer) [ alpha(real 0.05) ]
        . simulate reject = r(reject), reps(500) seed(12345): ///
        > syntax, n(24031) (.5)
        invalid syntax
        an error occurred when simulate executed syntax
        r(197);

        end of do-file

        Thanks for your reply,
        I have tried it and it comes up with the same error message

        Comment


        • #5
          Originally posted by daniel klein View Post
          I only read r(197) and the short program at the top of your post, so I won't comment on content. Syntax-wise, you want
          Code:
          syntax, n(integer) [ alpha(real 0.05) ]
          I was also making a mistake when running the program:
          Code:
          simulate reject = r(reject), reps(500) seed(12345): ///
          syntax, n(24031) (.5)
          Should have been:
          Code:
          simulate reject = r(reject), reps(500) seed(12345): ///
          sim_interaction_sacat, n(24031) alpha(.05)
          Its working now ! Thanks for fixing my mistake

          Comment

          Working...
          X