Help for Power analysis
"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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._