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
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

Comment