[R] Help for Power analysis

[Uwe Ligges]

"Ass.Prof. Nikom Thanomsieng" wrote:
> 
> Dear colleague,
>    I not sure this R code  is correctly ? 

Given the code you wrote is correct (didn't checked it, but looks OK in
the hurry):

> I would to show the number of Sample Size at  Sample Size Axis that 
> line draw from Power Axis (80%) from R code.
>    How I show this and select the most appropriate of
> this power (.79955687 - 80983575).

What does this (above) line mean?

I changed your code a little bit (loop and several lines are not
neccessary).
I put it into a function and added the lines to automatically draw the
lines it seems you want to be drawn (see below). Hope this helps...

> Thank for your help and answer.
> 
> Best Regards,
> Nikom Thanomsieng,
> Email: nikom at kku.ac.th
> 
> ....
> 
> #Power analysis: Sample size for Chi-Square 2x2, RxC: R software
> #Concept from SAS program to calculate power of ANOVA F-test
> #http://www2.tltc.ttu.edu/Westfall/images/5347/power_analysis_of_anova_f_test.htm
> 
> # and Cohen,J (1977). Statistical Power Analysis for Behavioral
> Sciences.
> #New York: Academic Press.
> #Asst. Prof. Nikom Thanomsieng. 29/09/2000
> #Department of Biostatistics & Demography. Faculty of Public Health.
> #Khon Kaen University. Thailand.
> #Email: nikom at kku.ac.th
> 
> #Modify data value of the first two line
> x1 <- c(6,9)
> x2 <- c(6,6)
> nc <-cbind(x1)
> nr  <-rbind(x2)
> data1 <- rbind(x1,x2)
> chi2<- chisq.test(data1,correct=F)$statistic
> Ntotal<-sum(data1)
> df<- ncol(nc-1)*nrow(nr-1)
> ifelse(df==1,w<- sqrt(chi2/Ntotal), w<-sqrt(chi2/(chi2+Ntotal)))
> Ntotal1<-900  #change this if power not enough
> alpha  <-0.05  #change this for One tailed =0.05
> ncp<-0
> chicrit<-NULL
> power<-NULL
> n<-NULL
> samplesize<-NULL
> for  (i in 1:Ntotal1){
>        ncp[i] <- w^2 * i
>        chicrit<-qchisq(1-alpha,df)
> power[i] <- 1-(pchisq(chicrit , df, ncp[i]))
> n[i]<-i
> samplesize<-cbind(n, ncp,power) }
> samplesize
> plot(n,power,type="l",col="red", lwd=1,
> panel.first = grid(10,10),
> main="Power as a function of Sample Size",
> xlab="Sample Size",
> ylab="Power" )
> segments(785, .8, 785,   0, col ="pink")
> segments(785, .8,    0, 0.8, col ="pink")
> mtext("Chi-Square", side = 3, line = 0.35,
>     outer = FALSE, at = mean(par("usr")[1:2]), cex = 1, font = 4)


power.analysis <- function(x1, x2, my.power = 0.8, 
    Ntotal1 = 900, alpha = 0.05){
  nc <- cbind(x1)
  nr <- rbind(x2)
  data1 <- rbind(x1, x2)
  chi2 <- chisq.test(data1, correct=FALSE)$statistic
  Ntotal <- sum(data1)
  df <- ncol(nc-1) * nrow(nr-1)
  w <- ifelse(df == 1, sqrt(chi2 / Ntotal), sqrt(chi2 / (chi2+Ntotal)))
  chicrit <- qchisq(1-alpha, df)
  n <- 1:Ntotal1
  ncp <- w^2 * n
  power <- 1 - (pchisq(chicrit , df, ncp))
  plot(n, power, type = "l", col = "red", lwd = 1, panel.first =
grid(10,10), 
    main = "Power as a function of Sample Size", xlab = "Sample Size", 
    ylab = "Power")
  min.samp <- which(power >= my.power)[1]
  segments(min.samp, my.power, min.samp,   0, col = "pink")
  segments(min.samp, my.power,   0, my.power, col = "pink")
  mtext("Chi-Square", side = 3, line = 0.35,
    outer = FALSE, at = mean(par("usr")[1:2]), cex = 1, font = 4)
  samplesize <- cbind(n, ncp, power)
  return(samplesize)
}

power.analysis(c(9, 6), c(6, 6))


Uwe Ligges
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._