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
r-help-bounces at r-project.org 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 at 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)
}