2

我写这段代码是为了得到 Dunnet anova post hoc test

import rpy2.robjects as ro
import rpy2.robjects.numpy2ri as npr
from rpy2.robjects.numpy2ri import numpy2ri as np2r
from rpy2.robjects.packages import importr
base = importr("base")
stats = importr('stats')
multcomp = importr('multcomp')

val = np2r(vFC)
exp = base.gl(4,6,24)
tiempo = base.factor(base.rep(base.c(0,2,5,10,15,30),4))

fmla = ro.Formula('val ~ tiempo + exp')
env = fmla.environment
env['val'] = val
env['tiempo'] = tiempo
env['exp'] = exp
anova = stats.aov(fmla)
print base.summary(anova)
Phoc = multcomp.glht(anova, linfct = ro.r('mcp(tiempo="Dunnet")'))
sPhoc = base.summary(Phoc)
print sPhoc

工作正常,但在分析后形成输出我得到源代码,我怎样才能摆脱该代码并仅从 Dunnet 分析中获得最终表。

            Df Sum Sq Mean Sq F value   Pr(>F)    
tiempo       5  91172   18234   8.788 0.000464 ***
exp          3  49402   16467   7.936 0.002108 ** 
Residuals   15  31125    2075                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: function (formula, data = NULL, projections = FALSE, qr = TRUE, 
    contrasts = NULL, ...) 
{
    Terms <- if (missing(data)) 
        terms(formula, "Error")
    else terms(formula, "Error", data = data)
    indError <- attr(Terms, "specials")$Error
    if (length(indError) > 1L) 
        stop(sprintf(ngettext(length(indError), "there are %d Error terms: only 1 is allowed", 
            "there are %d Error terms: only 1 is allowed"), length(indError)), 
            domain = NA)
    lmcall <- Call <- match.call()
    lmcall[[1L]] <- as.name("lm")
    lmcall$singular.ok <- TRUE
    if (projections) 
        qr <- lmcall$qr <- TRUE
    lmcall$projections <- NULL
    if (is.null(indError)) {
        fit <- eval(lmcall, parent.frame())
        if (projections) 
            fit$projections <- proj(fit)
        class(fit) <- if (inherits(fit, "mlm")) 
            c("maov", "aov", oldClass(fit))
        else c("aov", oldClass(fit))
        fit$call <- Call
        return(fit)
    }
    else {
        if (pmatch("weights", names(match.call()), 0L)) 
            stop("weights are not supported in a multistratum aov() fit")
        opcons <- options("contrasts")
        options(contrasts = c("contr.helmert", "contr.poly"))
        on.exit(options(opcons))
        allTerms <- Terms
        errorterm <- attr(Terms, "variables")[[1 + indError]]
        eTerm <- deparse(errorterm[[2L]], width.cutoff = 500L, 
            backtick = TRUE)
        intercept <- attr(Terms, "intercept")
        ecall <- lmcall
        ecall$formula <- as.formula(paste(deparse(formula[[2L]], 
            width.cutoff = 500L, backtick = TRUE), "~", eTerm, 
            if (!intercept) 
                "- 1"), env = environment(formula))
        ecall$method <- "qr"
        ecall$qr <- TRUE
        ecall$contrasts <- NULL
        er.fit <- eval(ecall, parent.frame())
        options(opcons)
        nmstrata <- attr(terms(er.fit), "term.labels")
        nmstrata <- sub("^`(.*)`$", "\\1", nmstrata)
        nmstrata <- c("(Intercept)", nmstrata)
        qr.e <- er.fit$qr
        rank.e <- er.fit$rank
        if (rank.e < length(er.fit$coefficients)) 
            warning("Error() model is singular")
        qty <- er.fit$residuals
        maov <- is.matrix(qty)
        asgn.e <- er.fit$assign[qr.e$pivot[1L:rank.e]]
        maxasgn <- length(nmstrata) - 1L
        nobs <- NROW(qty)
        if (nobs > rank.e) {
            result <- vector("list", maxasgn + 2L)
            asgn.e[(rank.e + 1):nobs] <- maxasgn + 1L
            nmstrata <- c(nmstrata, "Within")
        }
        else result <- vector("list", maxasgn + 1L)
        names(result) <- nmstrata
        lmcall$formula <- form <- update(formula, paste(". ~ .-", 
            deparse(errorterm, width.cutoff = 500L, backtick = TRUE)))
        Terms <- terms(form)
        lmcall$method <- "model.frame"
        mf <- eval(lmcall, parent.frame())
        xvars <- as.character(attr(Terms, "variables"))[-1L]
        if ((yvar <- attr(Terms, "response")) > 0L) 
            xvars <- xvars[-yvar]
        if (length(xvars)) {
            xlev <- lapply(mf[xvars], levels)
            xlev <- xlev[!sapply(xlev, is.null)]
        }
        else xlev <- NULL
        resp <- model.response(mf)
        qtx <- model.matrix(Terms, mf, contrasts)
        cons <- attr(qtx, "contrasts")
        dnx <- colnames(qtx)
        asgn.t <- attr(qtx, "assign")
        if (length(wts <- model.weights(mf))) {
            wts <- sqrt(wts)
            resp <- resp * wts
            qtx <- qtx * wts
        }
        qty <- as.matrix(qr.qty(qr.e, resp))
        if ((nc <- ncol(qty)) > 1) {
            dny <- colnames(resp)
            if (is.null(dny)) 
                dny <- paste0("Y", 1L:nc)
            dimnames(qty) <- list(seq(nrow(qty)), dny)
        }
        else dimnames(qty) <- list(seq(nrow(qty)), NULL)
        qtx <- qr.qty(qr.e, qtx)
        dimnames(qtx) <- list(seq(nrow(qtx)), dnx)
        for (i in seq_along(nmstrata)) {
            select <- asgn.e == (i - 1)
            ni <- sum(select)
            if (!ni) 
                next
            xi <- qtx[select, , drop = FALSE]
            cols <- colSums(xi^2) > 1e-05
            if (any(cols)) {
                xi <- xi[, cols, drop = FALSE]
                attr(xi, "assign") <- asgn.t[cols]
                fiti <- lm.fit(xi, qty[select, , drop = FALSE])
                fiti$terms <- Terms
            }
            else {
                y <- qty[select, , drop = FALSE]
                fiti <- list(coefficients = numeric(), residuals = y, 
                  fitted.values = 0 * y, weights = wts, rank = 0L, 
                  df.residual = NROW(y))
            }
            if (projections) 
                fiti$projections <- proj(fiti)
            class(fiti) <- c(if (maov) "maov", "aov", oldClass(er.fit))
            result[[i]] <- fiti
        }
        result <- result[!sapply(result, is.null)]
        class(result) <- c("aovlist", "listof")
        if (qr) 
            attr(result, "error.qr") <- qr.e
        attr(result, "call") <- Call
        if (length(wts)) 
            attr(result, "weights") <- wts
        attr(result, "terms") <- allTerms
        attr(result, "contrasts") <- cons
        attr(result, "xlevels") <- xlev
        result
    }
}(formula = val ~ tiempo + exp)

Linear Hypotheses:
            Estimate Std. Error t value Pr(>|t|)    
2 - 0 == 0     78.10      32.21   2.425  0.10346    
5 - 0 == 0    152.39      32.21   4.731  0.00113 ** 
10 - 0 == 0   140.06      32.21   4.348  0.00246 ** 
15 - 0 == 0   158.90      32.21   4.933  < 0.001 ***
30 - 0 == 0   180.73      32.21   5.611  < 0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported -- single-step method)
4

1 回答 1

3

如果我从高级接口使用字符串作为 R 代码,源代码就会消失

将此代码用于第三个块

ro.r.assign('val', val)
ro.r.assign('exp', exp)
ro.r.assign('tiempo', tiempo)
ro.r('anova <- aov(val ~ tiempo + exp)')
ro.r('s.anova <- summary(anova)')
ro.r('spHoc <- summary(glht(anova, linfct=mcp(tiempo="Dunnet")))')
print ro.r('s.anova')
print ro.r('spHoc')
ro.r('capture.output(s.anova, file = "anova.txt", append = TRUE)')
ro.r('capture.output(spHoc, file = "anova.txt", append = TRUE)')

返回没有源代码的结果

            Df Sum Sq Mean Sq F value   Pr(>F)    
tiempo       5  91172   18234   8.788 0.000464 ***
exp          3  49402   16467   7.936 0.002108 ** 
Residuals   15  31125    2075                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = val ~ tiempo + exp)

Linear Hypotheses:
            Estimate Std. Error t value Pr(>|t|)    
2 - 0 == 0     78.10      32.21   2.425  0.10306    
5 - 0 == 0    152.39      32.21   4.731  0.00110 ** 
10 - 0 == 0   140.06      32.21   4.348  0.00246 ** 
15 - 0 == 0   158.90      32.21   4.933  < 0.001 ***
30 - 0 == 0   180.73      32.21   5.611  < 0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported -- single-step method)

我仍然不知道为什么当我使用低级接口时 rpy2 的行为是这样的,但是现在它可以工作了

于 2012-06-28T23:01:23.867 回答