Hi below you can find a function I wrote to mimic a Minitab R+R procedure. It is intended for myself only so it is not commented. If you want to know how it functions contact me off-list.
Regards Petr [EMAIL PROTECTED] napsal dne 17.09.2007 05:55:00: > Dear List, > > Which is the most eficient method to perform a gage repeatability and > reproducibility using R? > > Thanks and best regards. > > Constant Depièreux > ______________________________________________ > R-help@r-project.org 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. #------------------------------------------------------------------------------------------------ # O&R analysis (Minitab) ## fit is result from aov model (vysledek~vzorek*operator) # operator is oprator number, vzorek is sample number ORanal<-function(operator,vzorek,vysledek, sigma=5.15, DM=NULL ,HM=NULL, toler=NULL, plotit=F, dig=4) { if (!is.factor(vzorek)) vzorek<-factor(vzorek) if (!is.factor(operator)) operator<-factor(operator) if (is.null(toler)) toler<-abs(HM-DM) fit<-aov(vysledek~vzorek*operator) volop<-nlevels(operator[,drop=T]) volvzor<-nlevels(vzorek[,drop=T]) opak<-length(operator)/(volop*volvzor) ctver<-anova(fit)[,3] result=NULL result[3]<-(ctver[3]-ctver[4])/opak result[5]<-ctver[4] result[1]<-(ctver[1]-ctver[4]-result[3]*opak)/(volop*opak) result[2]<-(ctver[2]-ctver[4]-result[3]*opak)/(volvzor*opak) if (result[2]<0) result[2]<-0 result[4]<-result[2]+result[3] result[6]<-result[4]+result[5] result[7]<-result[1]+result[6] result<-abs(result) proc<-result/result[7]*100 smodch<-sqrt(result) smodch5<-smodch*sigma proc.smodch<-smodch/smodch[7]*100 proc.toler<-smodch5/toler*100 result<-cbind(result,smodch,smodch5,proc,proc.smodch,proc.toler) result<-data.frame(result,row.names=c("Vzorky","Oper","Interakce","Reprod","Opak","O&R","Suma")) meno<-c("rozptyl","smer.odch", "5.15 sm.odch" ,"proc.rozpt." ,"proc.sm.odch", "toler.sm.odch") velik<-dim(result)[2] names(result)<-meno[1:velik] suma<-summary(fit) kateg<-trunc(smodch[1]/smodch[6]*1.41)+1 if (plotit) { par(mfrow=c(3,2)) posice<-volvzor*1:volop mat<-aggregate(vysledek,list(vzorek,operator),mean,na.rm=T) mat1<-aggregate(vysledek,list(vzorek,operator),sd,na.rm=T) #fig1 hmdm<-c(mean(mat$x)+(mean(mat1$x)/.9)/sqrt(opak)*3, mean(mat$x)-(mean(mat1$x)/.9)/sqrt(opak)*3) plot(mat$x,type="b",ylab="Prumery") abline(h=mean(mat$x),col=3) abline(h=hmdm,col=2) abline(v=posice+.5,col="grey",lwd=2,lty=2) text(posice-volvzor/2,max(mat$x),levels(operator[,drop=T])) #fig2 # mat<-reshape(mat,direction="wide",idvar="Group.1",timevar="Group.2") # matplot(as.matrix(mat[,2:(1+volop)]),type="b",pch=1:volop,col=1:volop,ylab="Prumer",axes=F) # axis(side=1,at=1:volvzor,labels=levels(vzorek[,drop=T])) # box() # axis(2) # legend(1,min(mat[,2:(1+volop)]),levels(operator[,drop=T]),pch=1:volop,col=1:volop,cex=.8,horiz=T,yjust=0,xjust=0) mat<-reshape(mat,direction="wide",idvar="Group.1",timevar="Group.2") matplot(as.matrix(mat[,2:(1+volop)]),type="b",pch=1:volop,col=1:volop,ylab="Prumer",axes=F) axis(side=1,at=1:volvzor,labels=levels(vzorek[,drop=T])) box() axis(2) legend(1,min(mat[,2:(1+volop)]),levels(operator[,drop=T]),pch=1:volop,col=1:volop,cex=.8,horiz=T,yjust=0,xjust=0) #fig3 hmdm<-c(mean(mat1$x)*sqrt(qchisq(.99865,opak,lower.tail=F)/(opak-1)), mean(mat1$x)*sqrt(qchisq(.99865,opak)/(opak-1))) plot(mat1$x,type="b",ylab="Smerod. odch") abline(h=mean(mat1$x),col=3) abline(h=hmdm,col=2) abline(v=posice+.5,col="grey",lwd=2,lty=2) text(posice-volvzor/2,max(mat1$x),levels(operator[,drop=T])) #fig4 boxplot(split(vysledek,operator),notch=T) #fig5 barplot(t(as.matrix(result[c(1,4,5,6),4:velik])),beside=T) #fig6 boxplot(split(vysledek,vzorek),notch=T) par(mfrow=c(1,1)) } rozlis<-smodch[6]*qnorm(0.975) list(tabulka=round(result,dig),vartab=suma,kategorie=kateg,rozlisitelnost=rozlis) } ______________________________________________ R-help@r-project.org 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.
