[R] Odp: GR&R - Best methods in R

[Petr PIKAL]

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)


}