[R] Code Help

[Stephen Beale]

I am trying to analyse some data and was given R code to do this with but
there seem to be errors in the code.  My level of knowledge is improving but
still limited.

The details are;

Data on clover lines; Lines.txt	attached.  Comma seperations

Code:

options(digits=3)

clover <- read.table("Lines.txt",header=T,sep=",")
vnames <- names(clover);nv <- length(vnames)
flevels <- levels(clover$Line)

par(oma=c(0,6,4,6),mfrow=c(2,1),mar=c(4,4,1,0))

clover$Y <- clover[,"LeafLmm"]
mn.tab <- tapply(clover$Y,list(clover$Line),mean)
sd.tab <- sqrt(tapply(clover$Y,list(clover$Line),var))

om <- order(mn.tab)
plot(mn.tab,sd.tab,las=1,xlab="mean",ylab="s.d.")
lines(lowess(mn.tab,sd.tab))


model.i <- glm(Y  ~ C(clover$Line,Cline),data=clover,family=gaussian),

Beta <- summary(model.i)$coefficients
new.order <- order(Beta[-1,1])  + 1
Obeta <- Beta[c(1,new.order),]

detect.diff <- quantile(Obeta,probs=seq(0.2,0.8,0.2))
std.dev <- (min(Obeta[-1,2]) + max(Obeta[-1,2]))/2*sqrt(5)
sink("ClovAttr.txt")
print(paste("_____________",vnames[(i+2)],"_____________"))
cat("__________ mean table _______ \n")
print(mn.tab[om])
cat("__________ s.d. table _______ \n")
print(sd.tab[om])
print(model.i$family)
print(Obeta)
print(detect.diff)
cat("_________________________________\n")
sink()

library(ctest)

stdev <- round(std.dev*c(1.25,1,0.75),2)
stderr <- round(stdev/sqrt(5),2)
samp.size <- 3:10
beta <- 0.8

ylab <- "detectable differences"
for (j in seq(along=stdev)){
   Delta <- numeric(0)
for(k in seq(along=samp.size)){
Power.calcs <- power.t.test(n=samp.size[k],delta=NULL,sd=stdev[j],
 sig.level=0.05,power=beta,type="one.sample",alternative="one.sided")

beta  <- Power.calcs$power
Delta <- c(Delta,Power.calcs$delta)
          }  #   end of the k loop
if(j==1){
max.Delta <- max(Delta)
plot(samp.size,Delta,type='l',xlab="n",ylab=ylab,las=1,
    ylim=c(0,max.Delta),xlim=c(3,10))
         }   #   end of if(j==1)
else lines(samp.size,Delta,lty=j)
         }   #   end of j loop
mtext(vnames[i],side=3,outer=T,cex=1.2)

legend(7,max.Delta,legend=paste("s.e",stderr),lty=1:length(stdev))
text(5,0.9*max.Delta,"power=0.8")

_________________________________
The issue is the number of leaves to be sampled.
As above I only needed to analyse as a linear model, Plot curves of
detectable difference @ 3 standard deviations etc.  All the code is there
but there seems to be errors.  model.i I don't know what Cline is supposed
to be (possibly clover$Line), the detect.diff line is incorrect, the looping
through J&K appears incorrect and I don't know what i in the vnames is
supposed to be.

Any help?  This is probably simple stuff for experienced people.

Stephen Beale

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: Lines.txt
Url: https://stat.ethz.ch/pipermail/r-help/attachments/20030521/20fd52eb/Lines.txt