R help - Mar 2007 - 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]]

Meesters, Erik wrote:
> 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 think your DF are upside-down:
> power.1<-1-pf(qf(.95,1,(T*R-2),ncp=0),1,(T*R-2),ncp=Dl)
> power.1
[1] 0.532651



-- 
   O__  ---- Peter Dalgaard             ?ster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907

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

Thierry.Onkelinx op inbo.be

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: r-help-bounces op stat.math.ethz.ch 
> [mailto:r-help-bounces op stat.math.ethz.ch] Namens Meesters, Erik
> Verzonden: woensdag 7 maart 2007 15:50
> Aan: r-help op stat.math.ethz.ch
> 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]]
> 
> ______________________________________________
> R-help op stat.math.ethz.ch 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.
>