Erik, I haven't seen an answer to your question, so I'll try to answer it. The problem is that you switched the degrees of freedom. You had:
1 - pf(qf(.95, Vl, 1, ncp = 0), Vl, 1, ncp = Dl) [1] 0.05472242 But it should be: 1 - pf(qf(.95, 1, Vl, ncp = 0), 1, Vl, ncp = Dl) [1] 0.532651 Cheers, Thierry ------------------------------------------------------------------------ ---- ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Reseach Institute for Nature and Forest Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology and quality assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 [EMAIL PROTECTED] www.inbo.be Do not put your faith in what statistics say until you have carefully considered what they do not say. ~William W. Watt A statistical analysis, properly conducted, is a delicate dissection of uncertainties, a surgery of suppositions. ~M.J.Moroney > -----Oorspronkelijk bericht----- > Van: [EMAIL PROTECTED] > [mailto:[EMAIL PROTECTED] Namens Meesters, Erik > Verzonden: woensdag 7 maart 2007 15:50 > Aan: [email protected] > Onderwerp: [R] Power calculation for detecting linear trend > > Dear people, > I've a problem in doing a power calculation. In Fryer and > Nicholson (1993), ICES J. mar. Sci. 50: 161-168 page 164 an > example is given with the following characteristics T=5, > points in time R=5, replicates > Var.within=0.1 > q=10, a 10% increase per year > The degrees of freedom for the test are calculated as > Vl=T*R-2=23 and the non-centrality parameter Dl=4.54. > Using this they get a power of 0.53, but the result that I'm > getting is 0.05472242. > > I've tried this several ways in R, but I'm not able to come > up with the same number. Am I doing something wrong in the > calculation of the power? > Here's my code: > > T<-5 > R<-5 > sigmasq<-0.1 > q<-10 > Vl<-(T*R)-2 > Dl<-(R*(T-1)*T*(T+1)/(12*sigmasq))*(log(1+(q/100)))^2 #Dl > result is still similar > > power.1<-1-pf(qf(.95,(T*R-2),1,ncp=0),(T*R-2),1,ncp=Dl) > > Thank you for any suggestions/help. > > I'm using R2.4.1, on windowsXP. > > Erik Meesters > > > > > [[alternative HTML version deleted]] > > ______________________________________________ > [email protected] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
