snoopyzhao

我现在遇到一个困难，肯请各位帮助，谢谢！



我有一个三因素的实验，采用方差分析（ANOVA）来分析各因素、及因素间交互作用的影响。在 R 中，一切都很简单。举一例来说明：

<br />
>      N <- c(0,1,0,1,1,1,0,0,0,1,1,0,1,1,0,0,1,0,1,0,1,1,0,0)<br />
>      P <- c(1,1,0,0,0,1,0,1,1,1,0,0,0,1,0,1,1,0,0,1,0,1,1,0)<br />
>      K <- c(1,0,0,1,0,1,1,0,0,1,0,1,0,1,1,0,0,0,1,1,1,0,1,0)<br />
>      yield <- c(49.5,62.8,46.8,57.0,59.8,58.5,55.5,56.0,62.8,55.8,69.5,<br />
+      55.0, 62.0,48.8,45.5,44.2,52.0,51.5,49.8,48.8,57.2,59.0,53.2,56.0)<br />
> <br />
>      npk <- data.frame(block = gl(6,4), N = factor(N), P = factor(P),<br />
+                        K = factor(K), yield = yield)<br />
>      npk.aov <- aov(yield ~ block + N + P + K + N:P + N:K + P:N, npk)<br />
>      npk.lm  <- lm(yield ~ block + N + P + K + N:P + N:K + P:N, npk)<br />
>      summary(npk.aov)<br />
            Df Sum Sq Mean Sq F value   Pr(>F)   <br />
block        5 343.29   68.66  4.8047 0.010458 * <br />
N            1 189.28  189.28 13.2459 0.002997 **<br />
P            1   8.40    8.40  0.5879 0.456914   <br />
K            1  95.20   95.20  6.6622 0.022808 * <br />
N:P          1  21.28   21.28  1.4893 0.244003   <br />
N:K          1  33.14   33.14  2.3188 0.151769   <br />
Residuals   13 185.77   14.29                    <br />
---<br />
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 <br />
>      summary(npk.lm)<br />
<br />
Call:<br />
lm(formula = yield ~ block + N + P + K + N:P + N:K + P:N, data = npk)<br />
<br />
Residuals:<br />
    Min      1Q  Median      3Q     Max <br />
-5.1583 -1.8250  0.1583  2.0000  4.3667 <br />
<br />
Coefficients:<br />
            Estimate Std. Error t value Pr(>|t|)    <br />
(Intercept)   51.683      2.559  20.195 3.36e-11 ***<br />
block2         3.425      2.673   1.281  0.22246    <br />
block3         6.750      2.673   2.525  0.02535 *  <br />
block4        -3.900      2.673  -1.459  0.16829    <br />
block5        -3.500      2.673  -1.309  0.21307    <br />
block6         2.325      2.673   0.870  0.40018    <br />
N1             9.850      2.673   3.685  0.00275 ** <br />
P1             0.700      2.183   0.321  0.75351    <br />
K1            -1.633      2.183  -0.748  0.46756    <br />
N1:P1         -3.767      3.087  -1.220  0.24400    <br />
N1:K1         -4.700      3.087  -1.523  0.15177    <br />
---<br />
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 <br />
<br />
Residual standard error: 3.78 on 13 degrees of freedom<br />
Multiple R-squared: 0.788,      Adjusted R-squared: 0.625 <br />
F-statistic: 4.833 on 10 and 13 DF,  p-value: 0.004942 <br />
<br />
> dummy.coef(npk.lm)<br />
Full coefficients are <br />
                                                                     <br />
(Intercept):     51.68333                                            <br />
block:                  1         2         3         4      5      6<br />
                    0.000     3.425     6.750    -3.900 -3.500  2.325<br />
N:                      0         1                                  <br />
                     0.00      9.85                                  <br />
P:                      0         1                                  <br />
                      0.0       0.7                                  <br />
K:                      0         1                                  <br />
                 0.000000 -1.633333                                  <br />
N:P:                  0:0       1:0       0:1       1:1              <br />
                 0.000000  0.000000  0.000000 -3.766667              <br />
N:K:                  0:0       1:0       0:1       1:1              <br />
                      0.0       0.0       0.0      -4.7              <br />
> <br />

我的问题是如何理解那个线性模型给出的系数，以及如何解释那个系数为什么是显著的？当 options(contrasts) 变化时，线性模型的系数会发生变化，如何理解这种变化？



谢谢！！

PCR, PLS 或者岭回归吧……

其实用方差分析没有问题啊。如果方差分析显示不同处理对观测值没有显著影响，那么就没有必要做多重比较了啊，也就是不存在用对照组与其余 9 组进行比较的前提了。如果有显著影响，这个时候可以根据需要采用 LSD, Duncan, TukeyHSD 等各种多重比较的程序，你只保留对照与其余的 9 组比较结果就可以，至于其它的两两比较，你可以不考虑嘛，呵呵……



至于方差不齐，按前面的说法，进行变换一般都可以满足……

Modern Applied Statistics With S



本版就有下载，呵呵……

到 MASS 中找找看看

install.packages("Hmisc")

看来也只能这样了，谢谢！

我没有见过 weibull 三参数概率分布是啥样子，不过如果有概率密度函数的表达式，那么一切都好办了，呵呵……

你是指 weibull 分布吗？我看了一下 r 的帮助，这个函数只受两个参数影响，即 shape 和 scale，你这里的 x, y, z 是指啥？



qw <- seq(0, 5, 0.01)

dw <- dweibull(qw, shape = 2, scale = 2)

plot(qw, dw, type = "l")

假定在 R 中已经画出一条曲线，那么如何画出这条曲线上某一指定点的切线？



非常感谢！

关于这一问题，可以看 R-help 上最近的讨论，似乎很热烈的样子……

用什么标准就选什么模型，没有什么优劣之分



另外，模型是用来解释问题的，你获得的模型可以解释你的问题，那就可以了，管它是什么呢……

还是看看 boot 包的文档吧，有现成的，为什么要自己写呢？

问题解决了，呵呵，我自己对问题的理解有问题

<br />
ahp.ri <- function(m,n)<br />
{<br />
    x <- m * (m - 1) * n<br />
    a <- c(1:9,1/(2:9))<br />
    y <- sample(a, x, replace = TRUE)<br />
    rec <- matrix(0, m, m)<br />
    ri <- numeric(length = n)<br />
    for (i in 1:n)<br />
    {<br />
        for (j in 1:m)<br />
        {<br />
            for (k in j:m)<br />
            {<br />
                if (j == k) <br />
                {<br />
                    rec[j,k] <- 1<br />
                }<br />
                else<br />
                {<br />
                    rec[j,k] <- y[1]<br />
                    y <- y[-1]<br />
                    rec[k,j] <- 1 / rec[j,k]<br />
                }<br />
            }<br />
        }<br />
        ri[i] <- eigen(rec)$values[1]<br />
    }<br />
    ri.ave <- (mean(ri) - m) / (m - 1) <br />
    return(ri.ave)<br />
}<br />

注意，不同的随机数产生方式，结果会有所不同

虽然 Saaty 的 AHP 书中给出了 平均随机一致性指标 RI 的值，也给出了它的计算方法（不是很详细，我看不到过程），根据我的理解，我写了一个 R 的函数，想重现他的结果，但结果却有很大差别。我想可能是我的理解有问题，但我也找不到他的方法的原始文献，这里的牛牛们给看看吧，谢谢了。

<br />
ahp.ri <- function(m,n)<br />
{<br />
    x <- m * (m - 1) * n<br />
    a <- c(1:9)<br />
    y <- sample(a, x, replace = TRUE)<br />
    rec <- matrix(0, m, m)<br />
    ri <- numeric(length = n)<br />
    for (i in 1:n)<br />
    {<br />
        for (j in 1:m)<br />
        {<br />
            for (k in 1:m)<br />
            {<br />
                if (j == k) rec[j,k] <- 1<br />
                if (j < k)<br />
                {<br />
                    rec[j,k] <- y[1]<br />
                    y <- y[-1]<br />
                }<br />
                if (j > k) rec[j,k] <- 1 / rec[k,j]<br />
            }<br />
        }<br />
        ri[i] <- (eigen(rec)$values[1] - m) / (m -1)<br />
    }<br />
    ri.ave <- mean(ri)<br />
    return(ri.ave)<br />
}<br />

运行结果，如，当阶数为 3，运行 500 次

<br />
ahp.ri(3,500)<br />
[1] 0.190372<br />

而 Saaty 的书中，这个值约为 0.58



请帮忙看看吧，谢谢啦！