Как получить R ^ 2 и P-значение для анализа R ELISA 4PL - PullRequest
Новая
       47

Как получить R ^ 2 и P-значение для анализа R ELISA 4PL

0 голосов
/ 08 января 2020

Я пытаюсь имитировать c анализ ELISA на графическом планшете, используя R, однако у меня возникают трудности с получением значения P и значения R ^ 2.

Я следовал учебному пособию: http://weightinginbayesianmodels.github.io/poctcalibration/calib_tut4_curve_ocon.html#unweighted -nonlinear-regression-in-r

Мне удалось получить большую часть необходимой информации, используя пакет с именем "minpack.lm", однако я не уверен, как подойти отсюда получаем значения R ^ 2 и P. 

  ODCalc1 <- c(.007, .072, .328, .988, 1.534, 1.983)
  ODCalc2 <- c(.006, .074, .361, .858, 1.612, 1.993)
  ODCalc <- (ODCalc1 + ODCalc2)/2

  concentration <- log10(c(1, 36, 180, 540, 1080, 1800))

  ocon <- data.frame(10^(concentration), "rep", ODCalc, stringsAsFactors = F)
  ocon$X.rep. <- as.numeric(ocon$X.rep.)
  ocon$X.rep. <- 1
  names(ocon) <- c("conc", "rep", "od")

  # Plot the O'Connell data
  par(mfrow = c(1, 2), cex.main = 1, mar = c(4, 4, 1, 2), oma = c(0.5, 0.5, 2.5, 0))
  plot(ocon$conc, ocon$od, pch = 21, bg = "grey", ylab = "Response (od)", 
       xlab = "Concentration")
  grid()
  # Plot on the log(x) scale
  plot(log(ocon$conc), ocon$od, pch = 21, bg = "grey", ylab = "Response (od)", 
       xlab = "log(concentration)")
  grid()
  title("O'Connell's ELISA: concentration on absolute (left) and log (right) scales",
        outer = T)

  par(mfrow = c(1, 1))

  # ------------ Function: 4PL curve function ---------------------------------  
  M.4pl <- function(x, small.x.asymp, inf.x.asymp, inflec, hill){
    f <- small.x.asymp + ((inf.x.asymp - small.x.asymp)/
                            (1 + (x / inflec)^hill))
    return(f)
  }
  # ------------- end ---------------------------------------------------------

  start.ocon <- c(small.x.asymp = 0.1, inf.x.asymp = 1, inflec = 3000, hill = -1)
  library(minpack.lm)
  uw.4pl <- nlsLM(od ~ M.4pl(conc, small.x.asymp, inf.x.asymp, inflec, hill), 
                  data = ocon,
                  start = start.ocon)
  data.4pl <- summary(uw.4pl)

  bottom.4pl <- data.4pl$parameters[1,1]
  top.4pl <- data.4pl$parameters[2,1]
  IC50.4pl <- data.4pl$parameters[3,1]
  HillSlope.4pl <- abs(data.4pl$parameters[4,1])

  RSS.p <- sum(residuals(uw.4pl)^2)
  TSS <- sum((ocon$od - mean(ocon$od))^2)
  r.squared <- 1-(RSS.p/TSS) # is this the proper way to get an r^2 value? It does not match what graphpad has which is an issue.

  # I have also read this should work, but since the model is a linear model instead of a Sigmoidal, 4PL, X is log (concentration) model
  model <- lm(concentration ~ poly(ODCalc, degree = 4, raw=T))
  summary(model) # R^2 is not the correct value I am looking for.

  # Not sure if sample data is needed but these were the values we were using to produce the values below
  sample.od.values1 <- c(0.275, 1.18, 0.085, 0.054, 0.119)
  sample.od.values2 <- c(0.263, 1.149, 0.068, 0.062, 0.109)
  sample.od.values <- (sample.od.values1+sample.od.values2)/2

Значения, подтверждающие методы, одинаковы:

bottom.4pl = 0.01657

top.4pl = 3.002

HillSlope = 1.222

R ^ 2 = 0.9978

R ^ 2 (скорректировано) = 0,9969

P-значение = 0,5106

Заранее благодарю за полезные советы!

1 Ответ

0 голосов
/ 09 января 2020

Так как R ^ 2 измеряет линейную ассоциацию, она обычно используется для линейной регрессии, но игнорируя, что это, кажется, дает нужные вам числа или, по крайней мере, числа, близкие к этим. Для скорректированной формулы R в квадрате см. https://en.wikipedia.org/wiki/Coefficient_of_determination#Adjusted_R2, а для значения p я предположил, что вы ищете значение p для гипотезы о том, что первый коэффициент равен нулю. 

RSS <- deviance(uw.4pl); RSS
## [1] 0.001514624

coef(uw.4pl) # coefficients/parameters
## small.x.asymp   inf.x.asymp        inflec          hill 
##    0.01654996    3.00261439 1033.53324214   -1.22171740 

R2 <- cor(ocon$od, fitted(uw.4pl))^2; R2
## [1] 0.9995529

n <- nobs(uw.4pl)    
p <- length(coef(uw.4pl))
adjR2 <- 1 - (1-R2) * (n - 1) / (n - p - 1); adjR2
## [1] 0.9977645

pvalue <- coef(summary(uw.4pl))[1, 4]; pvalue
## [1] 0.5486584
...