��Ū

���ۿޤ������������ʱߡʳ�Ψ�ʱߡˡ���ľ������ľ���ο���³�����MA��RMA �ˤ���ľ����������

����ˡ

scatter(x, y, ellipse=F, lrl=F, cb=F, ma=F, rma=F, alpha=0.05)

����

x        ��Ω�ѿ��ʲ�����
y        ��°�ѿ��ʽļ���
ellipse  ��Ψ�ʱߤ������Ȥ� ellipse=TRUE ����� ��
lrl      ��ľ���������Ȥ� lel=TRUE ����� ��
cb       ��ľ���ο���³��Ӥ������Ȥ� cb=TRUE �����
ma       Major Axis regression �ˤ���ľ���������Ȥ� ma=TRUE ����� ��
rma      Reduced Major Axis regression �ˤ���ľ���������Ȥ� rma=TRUE ����� ��
alpha    ��Ψ�ʱߡ���ľ���ο���³��Ӥ������Ȥ��Φ��������١�(1-���ˡ�100%
acc      ��Ψ�ʱߤΤʤ�餫��
xlab     x ��̾��
ylab     y ��̾��

����TRUE/FALSE �ʳ��ˡ�Ĺ�� 3 �ο��ͥ٥��ȥ�ǡ�
 (����, ����, ��) �����Ǥ��롣���ͤ� lty, lwd, col �򻲾ȤΤ���
�����ʱߤΤȤ��ϡ����ϡ�border �ο�

������

���󥹥ȡ���ϡ��ʲ��� 1 �Ԥ򥳥ԡ�����R ���󥽡���˥ڡ����Ȥ���
source("http://aoki2.si.gunma-u.ac.jp/R/src/scatter.R", encoding="euc-jp")

# ���ۿޤ������������ʱߡʳ�Ψ�ʱߡˡ���ľ������ľ���ο���³��ӡ�MA��RMA �ˤ���ľ��������
scatter <- function( x,                                      # ��Ω�ѿ��ʲ����ˤȤ��
                        y,                                      # ��°�ѿ��ʽļ��ˤȤ��
                        ellipse=FALSE,                          # ��Ψ�ʱߤ�ä���Ȥ��� TRUE *1
                        lrl=FALSE,                              # ��ľ����ä���Ȥ��� TRUE *1
                        cb=FALSE,                               # ��ľ���ο���³���ä���Ȥ��� TRUE
                        ma=FALSE,                               # MA �ˤ���ľ����ä���Ȥ��� TRUE *1
                        rma=FALSE,                              # RMA �ˤ���ľ����ä���Ȥ��� TRUE *1
                        alpha=0.05,                             # 1-�� �������١ʿ���Ψ��
                        acc=2000,                               # ��Ψ�ʱߤΤʤ�餫��
                        xlab=NULL,                              # x ��̾��
                        ylab=NULL)                              # y ��̾��
                                                                # *1 TRUE/FALSE �ʳ��ˡ�Ĺ�� 3 �ο��ͥ٥��ȥ�ǡ�
                                                                # (����, ����, ��) �����Ǥ��롣���ͤ� lty, lwd, col �򻲾ȤΤ���
{                                                               # �����ʱߤΤȤ��ϡ����ϡ�border �ο�
        comp <- function(a)
        {
                if (length(a) == 1) a[2] <- a[3] <- 1
                if (length(a) == 2) a[3] <- 1
                return(a)
        }       
        
        ellipse.draw <- function(x, y, alpha, acc, ellipse,  # �����ʱߤ�����
                                 xlab, ylab)
        {
                ellipse <- comp(ellipse)
                vx <- var(x)
                vy <- var(y)
                vxy <- var(x, y)
                lambda <- eigen(var(cbind(x, y)))$values
                a <- sqrt(vxy^2/((lambda[2]-vx)^2+vxy^2))
                b <- (lambda[2]-vx)*a/vxy
                theta <- atan(a/b)
                k <- sqrt(-2*log(alpha))
                l1 <- sqrt(lambda[1])*k
                l2 <- sqrt(lambda[2])*k
                x2 <- seq(-l1, l1, length.out=acc)
                tmp <- 1-x2^2/l1^2
                y2 <- l2*sqrt(ifelse(tmp < 0, 0, tmp))
                x2 <- c(x2, rev(x2))
                y2 <- c(y2, -rev(y2))
                s0 <- sin(theta)
                c0 <- cos(theta)
                xx <- c0*x2+s0*y2+mean(x)
                yy <- -s0*x2+c0*y2+mean(y)
                rngx <- range(c(x, xx))
                rngy <- range(c(y, yy))
                plot(rngx, rngy, type="n", xlab=xlab, ylab=ylab)
                polygon(xx, yy, lty=ellipse[1], lwd=ellipse[2], border=ellipse[3])
        }

        conf.limit <- function(x, y, cb, alpha, lrl)         # ��ľ���ȿ���³��Ӥ�����
        {
                lrl <- comp(lrl)
                n <- length(x)
                b <- var(x, y)/var(x)
                a <- mean(y)-b*mean(x)
                abline(a, b, lty=lrl[1], lwd=lrl[2], col=lrl[3])
                if (cb) {
                        sx2 <- var(x)*(n-1)
                        R <- par()$usr                               # �������ϰ�
                        x1 <- seq(R[1], R[2], length.out=2000)
                        y1 <- a+b*x1
                        ta <- -qt(alpha/2, n-2)
                        Ve <- (var(y)-var(x, y)^2/var(x))*(n-1)/(n-2)
                        temp <- ta*sqrt(Ve)*sqrt(1/n+(x1-mean(x))^2/sx2)
                        y2 <- y1-temp
                        lines(x1, y2, lty="dotted")
                        y2 <- y1+temp
                        lines(x1, y2, lty="dotted")
                        temp <- ta*sqrt(Ve)*sqrt(1+1/n+(x1-mean(x))^2/sx2)
                        y2 <- y1-temp
                        lines(x1, y2, lty="dashed")
                        y2 <- y1+temp
                        lines(x1, y2, lty="dashed")
                }
                cat("LS(least squares)--------------------\n")
                list(intercept=a, slope=b)
        }
        
        MA <- function(x, y, ma)                             # Major Axis regression
        {
                ma <- comp(ma)
                s2 <- cov(cbind(x, y))
                b <- s2[1, 2]/(eigen(s2)$values[1]-s2[2, 2])
                a <- mean(y)-b*mean(x)
                abline(a, b, lty=ma[1], lwd=ma[2], col=ma[3])
                cat("MA(major axis)--------------------\n")
                list(intercept=a, slope=b)
        }
        
        RMA <- function(x, y, rma)                                   # Reduced Major Axis regression
        {
                rma <- comp(rma)
                b <- sign(cor(x, y))*sqrt(var(y)/var(x))
                a <- mean(y)-b*mean(x)
                abline(a, b, lty=rma[1], lwd=rma[2], col=rma[3])
                cat("RMA(reduced major axis)--------------------\n")
                list(intercept=a, slope=b)
        }

        if (is.null(xlab)) xlab <- deparse(substitute(x))    # x ��̾��
        if (is.null(ylab)) ylab <- deparse(substitute(y))    # y ��̾��
        OK <- complete.cases(x, y)                           # ��»�ͤ���ĥ����������
        x <- x[OK]
        y <- y[OK]
        if (ellipse[1]) {                                       # �ǡ����ݥ���Ȥ�ޡ������ƴ����ʱߤ�����
                ellipse.draw(x, y, alpha, acc, ellipse, xlab, ylab)
                points(x, y)
        }
        else {                                                  # ���ۿޤΤ�
                plot(x, y, xlab=xlab, ylab=ylab)
        }
        if (lrl[1]) {                                           # ��ľ���ȿ���³���
                print(conf.limit(x, y, cb, alpha, lrl))
        }
        if (ma[1]) {                                            # Major Axis regression
                print(MA(x, y, ma))
        }
        if (rma[1]) {                                           # Reduced Major Axis regression
                print(RMA(x, y, rma))
        }
}


������

> x <- c(132, 146, 140, 196, 132, 154, 154, 168, 140, 140, 156, 114, 134, 116, 150, 178, 150, 120, 150, 146)
> y <- c(90, 90, 84, 96, 90, 90, 74, 92, 60, 82, 80, 62, 80, 80, 76, 98, 86, 70, 80, 80)

> scatter(x, y, el=TRUE, lrl=TRUE, cb=TRUE, ma=TRUE, rma=TRUE)	# �������������������Ƥߤ�
LS(least squares)--------------------	# ľ���󵢤����Ҥȷ���
$intercept
[1] 36.90722

$slope
[1] 0.3092784

MA(major axis)--------------------	# Major Axis regression �ˤ�����Ҥȷ���
$intercept
[1] 28.95068

$slope
[1] 0.3638499

RMA(reduced major axis)--------------------	# Reduced Major Axis regression �ˤ�����Ҥȷ���
$intercept
[1] 7.29765

$slope
[1] 0.5123618

�ʲ��Τ褦�ʿޤ��������ʷ����εޤʤ�Τ����ˡ�RMA, MA, ���̤β�ľ����



> scatter(x, y, ellipse=c(1,1,6), lrl=c(1,2,1), ma=c(2,2,2), rma=c(3,3,4)) # ������������λ����Ȥǻ��ꤹ��Ȥ�



 ����ڡ���