English 中文(简体)
从适合R的模型中获取未经调整的价值
原标题:Getting unadjusted values from model fit in R
  • 时间:2011-11-23 09:45:06
  •  标签:
  • r

我想有一个经过调整和未经调整(精减)的表格。 这似乎是一个共同的问题,但我可以指出,如何这样做,而不要再重复每个变量的倒退。

这里我如何获得调整的数值:

library("survival")
library("timereg")
data(sTRACE)

# Basic cox regression    
surv <- with(sTRACE, Surv(time/365,status==9))
fit1 <- coxph(surv~age+sex+diabetes+chf+vf, data=sTRACE)
x <- cbind(exp(coef(fit1)), exp(confint(fit1)))
test <- apply(x, 1, FUN=function(x){
    x <- round(x, 1)
    txt <- paste(x[1], " (95% CI ", min(x[2:3]), "-", max(x[2:3]), ")", sep="")
    return(txt)
})
test

测试变量现在是病媒:

> test
                   age                    sex               diabetes 
"1.1 (95% CI 1.1-1.1)" "1.4 (95% CI 1.1-1.9)"   "1.5 (95% CI 1-2.2)" 
                   chf                     vf 
"2.1 (95% CI 1.6-2.8)" "2.3 (95% CI 1.4-3.8)" 

I woul like to add this to a 2-dimensional matrix where one column is the unadjusted value something like this:

                       Adjusted             Unadjusted
     age "1.1 (95% CI 1.1-1.1)" "1.1 (95% CI 1.1-1.1)"

2. 如果产生不调整的价值,则:

fit2 <- coxph(surv~age, data=sTRACE)
x <- cbind(exp(coef(fit2)), exp(confint(fit2)))
test <- apply(x, 1, FUN=function(x){
    x <- round(x, 1)
    txt <- paste(x[1], " (95% CI ", min(x[2:3]), "-", max(x[2:3]), ")", sep="")
    return(txt)
})
test

这可能是通过更新职能来实现的,但我认为应当有一些自动化过程,因为这是一个共同的做法。


UPDATE

After some thinking and with some inspiration from the answer I wrote this function:

print_adjusted_and_unadjusted <- function(fit, digits=2){
    # Just a prettifier for the output an alternative could be:
    # paste(round(x[,1],1), " (95% CI ", min(round(x[,2:3])), "-", max(round(x[,2:3])), ")", sep="")  
    get_coef_and_ci <- function(fit){
        # Just to make sure that it gives 1.0 and 
        # not 1 if digits = 1, in cases where a 
        # adding another decimal that is used
        # since everyone is so hyped about p-val < 0.05
        add_zero_to_var <- function(x){
            ret <- round(as.double(x), digits)
            if (x == 1){
                ret <- round(x, digits+1)
                if (ret == 1){
                    ret <- paste("1.", paste(rep("0", digits), collapse=""), sep="")
                }
            }else if(nchar(as.character(x)) < digits + 2){
                add_zeros <- digits + 2 - nchar(as.character(x))
                ret <- paste(x, paste(rep("0", add_zeros), collapse=""), sep="")
            }
            return(ret)
        }

        # Get coefficients and conf. interval
        my_coefficients <- coef(fit)
        ci <- confint(fit)

        # Use the exp() if logit or cox regression
        if ("coxph" %in% class(fit) ||
            ("glm" %in% class(fit) &&
            fit$family$link == "logit")){
            my_coefficients <- exp(my_coefficients)
            ci <- exp(ci)
        }

        if (length(my_coefficients) > 1){
            my_coefficients <- tapply(my_coefficients, 1:length(my_coefficients), FUN = add_zero_to_var)
        }else{
            my_coefficients <- add_zero_to_var(my_coefficients)
        }

        ci <- apply(ci, 1, FUN=function(x){
                ci <- round(x, digits)
                for(i in 1:2){
                    ci[i] <- add_zero_to_var(ci[i])
                }
                return(paste(ci[1], "-", ci[2], sep=""))
            })
        ret_val <- cbind(my_coefficients, ci)
        colnames(ret_val) <- c("", "2.5% - 97.5%")
        rownames(ret_val) <- names(coef(fit))
        return(ret_val)
    }

    # Extract all the term names
    all.terms <- terms(fit) 
    var_names <- attr(all.terms,  term.labels )

    # Skip variables consisting of
    # functions such as spline, strata variables
    regex_for_unwanted_vars <- "^(strat[a]{0,1}|ns|rcs|bs|pspline)[(]"
    skip_variables <- grep(regex_for_unwanted_vars, var_names)

    # Get the adjusted variables
    adjusted <- get_coef_and_ci(fit)
    # When using splines, rcs in cox regression this shows a little different

    # Remove all the splines, rcs etc
    rn <- rownames(adjusted)
    remove_1 <- grep("( {1,}|[[][0-9]+[]]|[)][0-9]+)$", rn)
    remove_2 <- grep("^(strat[a]{0,1}|ns|rcs|bs)[(]", rn)
    adjusted <- adjusted[-union(remove_1, remove_2), ]
    if ("cph" %in% class(fit)){
        remove_3 <- grep("^rcs[(]", var_names)
        adjusted <- adjusted[-remove_3, ]
    }

    unadjusted <- c() 
    for(variable in var_names[-skip_variables]){
        interaction_variable <- length(grep(":", variable)) > 0

        # If it s an interaction variable the
        # interacting parts have to be included  
        if (interaction_variable){
            variable <- paste(paste(unlist(strsplit(variable, ":")), sep="", collapse=" + "), variable, sep=" + ")
        }

        # Run the same fit but with only one variable
        fit_only1 <- update(fit, paste(".~", variable))

        # Get the coefficients processed with some advanced 
        # round part()
        new_vars <- get_coef_and_ci(fit_only1)

        # If interaction then we should only keep the 
        # interaction part - the other variables are
        # always included by default and need therefore
        # to be removed
        if (interaction_variable){
            new_vars <- new_vars[grep("[*:]", rownames(new_vars)),]
        }

        # Add them to the previous
        unadjusted <- rbind(unadjusted, new_vars)
    }

    # If regression contains (Intercept) 
    # this is meaningless for the comparison
    # of adjusted and unadjusted values 
    if ("(Intercept)" %in% rownames(unadjusted)){
        unadjusted <- unadjusted[rownames(unadjusted) != "(Intercept)", ]
        unadjusted <- rbind(c("-", "-"), unadjusted)
        rownames(unadjusted)[1] <- "(Intercept)"
    }

    both <- cbind(unadjusted, adjusted)
    colnames(both) <- c("Unadjusted", "95% CI", "Adjusted", "95% CI")
    return(both)
}

它给我一个四维阵列:

    Unadjusted 95% CI      Adjusted 95% CI     
age "0.74"     "0.68-0.81" "0.69"   "0.62-0.76"
....

我将此事与表(或晚间)一并使用。

xtable(print_adjusted_and_unadjusted(fit.oa.base.model), align="lrcrc")

我在10(a)、c(b)和cxph(b)上测试了它,似乎还在工作。

感谢你们的帮助,并希望这部法律不仅仅供我使用。

最佳回答

首先,鉴于你从<条码>cbind上取得了成绩,你不需要<条码>> 应用程序<>>>。

test<-paste(round(x[,1],1), " (95% CI ", min(round(x[,2:3])), "-", max(round(x[,2:3])), ")", sep="")

这应当产生同样的结果。

现在,如果你想要掌握一些不同的变数,那么你就不得不把你的公式作为特性(注:我假定测试是这里完全模式的结果,按照你的代码,这样我就可以使用姓名):

unadjusted<-sapply(names(test), function(curname){
  curfrm<-paste("surv", curname, sep="~")
  curfit<-coxph(as.formula(curfrm), data=sTRACE)
  curx <- cbind(exp(coef(fit1)), exp(confint(fit1)))
  paste(round(curx[,1],1), " (95% CI ", min(round(curx[,2:3])), "-", max(round(curx[,2:3])), ")", sep="")
})

Now you can just cbind test en unadjusted for your desired effect.

问题回答

暂无回答




相关问题
How to plot fitted model over observed time series

This is a really really simple question to which I seem to be entirely unable to get a solution. I would like to do a scatter plot of an observed time series in R, and over this I want to plot the ...

REvolution for R

since the latest Ubuntu release (karmic koala), I noticed that the internal R package advertises on start-up the REvolution package. It seems to be a library collection for high-performance matrix ...

R - capturing elements of R output into text files

I am trying to run an analysis by invoking R through the command line as follows: R --no-save < SampleProgram.R > SampleProgram.opt For example, consider the simple R program below: mydata =...

R statistical package: wrapping GOFrame objects

I m trying to generate GOFrame objects to generate a gene ontology mapping in R for unsupported organisms (see http://www.bioconductor.org/packages/release/bioc/vignettes/GOstats/inst/doc/...

Changing the order of dodged bars in ggplot2 barplot

I have a dataframe df.all and I m plotting it in a bar plot with ggplot2 using the code below. I d like to make it so that the order of the dodged bars is flipped. That is, so that the bars labeled "...

Strange error when using sparse matrices and glmnet

I m getting a weird error when training a glmnet regression. invalid class "dgCMatrix" object: length(Dimnames[[2]]) must match Dim[2] It only happens occasionally, and perhaps only under larger ...

Generating non-duplicate combination pairs in R

Sorry for the non-descriptive title but I don t know whether there s a word for what I m trying to achieve. Let s assume that I have a list of names of different classes like c( 1 , 2 , 3 , 4 ) ...

Per panel smoothing in ggplot2

I m plotting a group of curves, using facet in ggplot2. I d like to have a smoother applied to plots where there are enough points to smooth, but not on plots with very few points. In particular I d ...

热门标签