如何在r中使用Predict.lm来反转回归-Java 学习之路

我在数据帧calvarbyruno.1中有一些数据，其中变量Nominal和PAR表示使用特定分析技术分析一组标准时发现的峰面积比（PAR），以及该数据的两个lm模型（线性和二次）对于PAR~Nominal的关系 . 我正在尝试使用predict.lm函数来反馈计算Nominal值，给定我的PAR值，但是既可以预测.lm和fit也只能给出PAR值 . 我慢慢失去了我的魔力，任何人都可以帮忙吗？

calvarbyruno.1数据帧

structure(list(Nominal = c(1, 3, 6, 10, 30, 50, 150, 250), Run = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("1", "2", "3"), class = "factor"), 
    PAR = c(1.25000000000000e-05, 0.000960333333333333, 0.00205833333333334, 
    0.00423333333333333, 0.0322333333333334, 0.614433333333334, 
    1.24333333333333, 1.86333333333333), PredLin = c(-0.0119152187070942, 
    0.00375925114245899, 0.0272709559167888, 0.0586198956158952, 
    0.215364594111427, 0.372109292606959, 1.15583278508462, 1.93955627756228
    ), PredQuad = c(-0.0615895732702735, -0.0501563307416599, 
    -0.0330831368244257, -0.0104619953693943, 0.100190275883806, 
    0.20675348710041, 0.6782336426345, 1.04748729725370)), .Names = c("Nominal", 
"Run", "PAR", "PredLin", "PredQuad"), row.names = c(NA, 8L), class = "data.frame")

线性模型

summary(callin.1)

Call:
lm(formula = PAR ~ Nominal, data = calvarbyruno.1, weights = Nominal^calweight)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0041172 -0.0037785 -0.0003605  0.0024465  0.0071815 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.007083   0.005037  -1.406   0.2093  
Nominal      0.005249   0.001910   2.748   0.0334 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.004517 on 6 degrees of freedom
Multiple R-squared: 0.5572,     Adjusted R-squared: 0.4835 
F-statistic: 7.551 on 1 and 6 DF,  p-value: 0.03338

二次模型

> summary(calquad.1)

Call:
lm(formula = PAR ~ Nominal + I(Nominal^2), data = calvarbyruno.1)

Residuals:
        1         2         3         4         5         6         7         8 
 0.053366  0.033186  0.002766 -0.036756 -0.211640  0.177012 -0.021801  0.003867 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -6.395e-02  6.578e-02  -0.972  0.37560   
Nominal       1.061e-02  2.205e-03   4.812  0.00483 **
I(Nominal^2) -1.167e-05  9.000e-06  -1.297  0.25138   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.128 on 5 degrees of freedom
Multiple R-squared: 0.9774,     Adjusted R-squared: 0.9684 
F-statistic: 108.2 on 2 and 5 DF,  p-value: 7.658e-05

但Predict给了我这些 Value ，这两个看起来都是错的（虽然我无法弄清楚它在做什么，这对第二套不同？

> predict(callin.1)
           1            2            3            4            5            6 
-0.001834123  0.008663451  0.024409812  0.045404959  0.150380698  0.255356437 
           7            8 
 0.780235132  1.305113826 
> predict(callin.1,type="terms")
      Nominal
1 -0.32280040
2 -0.31230282
3 -0.29655646
4 -0.27556131
5 -0.17058558
6 -0.06560984
7  0.45926886
8  0.98414755
attr(,"constant")
[1] 0.3209663

编辑：正如已经指出的那样，我对自己要实现的目标并不十分清楚，所以我会尝试更好地自我解读 .

该数据来自一组已知浓度标准（标称值）的分析，该标准给出了一组特定的响应或峰面积比（PAR） . 我想展示哪种模型最适合这些数据，然后分析未知样品以找出它们的浓度 .

我正在努力跟随为此工作的其他人，这涉及到;
a）通过找到PAR的内部运行方差并将其拟合到log的模型（方差（PAR））=博客（标称值），找到要使用的适当权重，其中B将是要使用的权重（舍入到最近）整数）
b）将每次运行的数据拟合为线性模型（PAR = a bNominal）和二次模型（PAR = A Bnominal cNominal ^ 2）
c）返回计算每个标准的发现浓度，并与标称浓度进行比较，得出偏差
d）评估校准范围内的偏差，并根据偏差选择模型

这个问题试图做c） . R邮件列表的帖子表明，仅使用反转术语进行回归是不合适的，我可以手动对线性模型进行计算，但我正在与二次模型进行斗争 . 似乎从搜索R邮件列表中看到其他人想要做同样的事情 .

1 回答

好吧，我实际上不得不尝试这个，在看了各种各样的事情之后我写了一个函数来找到二次方程的根 .

invquad<-function(a,b,c,y,roots="both", xmin=(-Inf), xmax=(Inf),na.rm=FALSE){
#Calculate the inverse of a quadratic function y=ax^2+bx+c (ie find x when given y)
#Gives NaN with non real solutions
root1<-sqrt((y-(c-b^2/(4*a)))/a)-(b/(2*a))
root2<--sqrt((y-(c-b^2/(4*a)))/a)-(b/(2*a))
if (roots=="both") {
    root1<-ifelse(root1<xmin,NA,root1)  
    root1<-ifelse(root1>xmax,NA,root1)  
    root2<-ifelse(root2<xmin,NA,root2)  
    root2<-ifelse(root2>xmax,NA,root2)      
    result<-c(root1,root2)
    if (na.rm) result<-ifelse(is.na(root1),root2, result)
    if (na.rm) result<-ifelse(is.na(root2),root1,result)
    if (na.rm) result<-ifelse(is.na(root1)&is.na(root2),NA,result)
},roots="both"
if (roots=="min")
    result<-pmin(root1,root2, NA.rm=TRUE)
if (roots=="max")
    result<-pmax(root1,root2, NA.rm=TRUE)
result
}

所以，鉴于原始数据

> PAR
[1] 0.0000125000 0.0009603333 0.0020583333 0.0042333333 0.0322333333 0.6144333333
[7] 1.2433333333 1.8633333333
> Nominal
[1]   1   3   6  10  30  50 150 250

我们可以进行分析，找到系数然后找到反函数，对我们期望的标称值使用一些合理的限制...

lm(PAR~Nominal+I(Nominal^2))->bob
> bob[[1]][[3]]
[1] -1.166904e-05 # Nominal^2
> bob[[1]][[2]]
[1] 0.01061094 # Nominal
> bob[[1]][[1]]
[1] -0.06395298 # Intercept
> invquad(bob[[1]][[3]],bob[[1]][[2]],bob[[1]][[1]],y=PAR,xmin=-0.2,xmax=300,na.rm=TRUE)
[1]   6.068762   6.159306   6.264217   6.472106   9.157041  69.198703 146.949154
[8] 250.811211

希望这可以帮助....

回复于 2024-05-05T15:04:50+08:00

如何在r中使用Predict.lm来反转回归

1 回答

相关问题