Announcement

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

  • powersim

    I am having trouble getting the simulations to work for a particular problem and was wondering if anybody has already looked into this?

    The problem is that we have got a sample size based on a difference in means. We are going to use the negative binomial distribution to look at the difference between the number of events over a period of 2 years. One of the independent committee members thinks that the power we will have will be much smaller than we think (based on the difference in means)

    I was wondering if anybody has done something similar with powersim before and if so could you give me a pointer?

  • #2
    Have you looked at:
    https://www.stata.com/support/faqs/s...by-simulation/

    Other than that we have no information to tell you what is wrong with the program you did not show us.
    ---------------------------------
    Maarten L. Buis
    University of Konstanz
    Department of history and sociology
    box 40
    78457 Konstanz
    Germany
    http://www.maartenbuis.nl
    ---------------------------------

    Comment


    • #3
      This is the program I have so far. But it doesn't work.
      tony


      forvalues i = 0/1 {

      if `i' == 1 local robust robust


      powersim , ///
      b(0.33) ///
      pos(1) ///
      sample(200(100)600) ///
      nreps(4000) ///
      family(gaussian) ///
      link(identity) ///
      cov1(x1 _bp block 2) ///
      saving(psim_data`i', replace) ///
      dofile(simulate_dofile, replace) ///
      : glm y x1, family(nbinomial) link(log) `robust'
      }
      gen method = "ROBUST"
      append using psim_data0.dta, nolabel
      quietly replace method = "OIM" if method == ""
      }

      Comment


      • #4
        I think you are better off doing this from first principle

        Code:
        // effect is in terms of incidence rate ratio not in terms of raw coefficients
        clear all
        
        program define sim, rclass
            syntax, delta(real) effect(real) n(real) [robust]
            
            drop _all
            set obs `n'
            gen x = _n <= _N/2
            gen mu = exp(.5+ln(`effect')*x)
            gen lambda = rgamma(`delta' , mu)
            gen y = rpoisson(lambda)
            nbreg y x, irr `robust'
            return scalar p = 2*normal(-abs(_b[x] /_se[x]))
        end
        
        // loop over different scenarios
        tempname memhold
        tempfile results
        postfile `memhold' n robust power using `results'
        forvalues n = 200(100)600{
            forvalues rob = 0/1 {
                if `rob' == 1 {
                    local robust = "robust"
                }
                else {
                    local robust = ""
                }
                simulate p = r(p), reps(4000) nodots : sim, delta(2) effect(1.25) n(`n') `robust'
                count if p < .05
                post `memhold' (`n') (`rob') (r(N)/_N)
            }
        }
        postclose `memhold'
        
        // display results
        use `results', clear
        label define robust 0 "OIM" 1 "robust"
        label value robust robust
        list, sepby(n) noobs
        ---------------------------------
        Maarten L. Buis
        University of Konstanz
        Department of history and sociology
        box 40
        78457 Konstanz
        Germany
        http://www.maartenbuis.nl
        ---------------------------------

        Comment


        • #5
          thank you so much for your help
          I will give that a go and let you know how I get on.
          thanks again for this you are a star!!

          Comment

          Working...
          X