"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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._