Comment utiliser Predict.lm dans r pour inverser la régression

Question

J'ai des données dans un cadre de données calvarbyruno.1 avec des variables Nominal et PAR qui représentent le rapport de surface de pic (PAR) trouvé à partir de l'analyse d'un ensemble de normes utilisant une technique analytique particulière, et deux modèles lm de ces données (linéaire et quadratique) pour la relation PAR ~ Nominal. J'essaie d'utiliser la fonction predict.lm pour calculer à rebours les valeurs nominales, compte tenu de mes valeurs PAR, mais predict.lm et fitted semblent tous deux ne me donner que des valeurs PAR. Je suis en train de perdre mon mojo, quelqu'un peut-il m'aider ?

calvarbyruno.1 cadre de données

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")

Modèle linéaire

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 

Modèle quadratique

> 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

Mais Predict me donne ces valeurs, qui semblent toutes deux erronées (bien que je n'arrive pas à comprendre ce qu'il fait de différent pour le deuxième ensemble ?

> 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

EDIT : Comme cela a été souligné, je n'ai pas été très clair sur ce que j'essaie d'obtenir, je vais donc essayer de mieux m'exprimer.

Les données proviennent de l'analyse d'un ensemble d'étalons de concentrations connues (nominales) qui donne un ensemble particulier de réponses, ou rapports de surface de pic (PAR). Je veux montrer quel modèle correspond le mieux à ces données et l'utiliser pour analyser des échantillons inconnus afin de déterminer leur concentration.

J'essaie de suivre quelqu'un d'autre qui travaille pour ça, ce qui implique ;
a) Trouver le poids approprié à utiliser, en trouvant la variance intra-série de la RAP et en l'adaptant à un modèle de log(Variance(RAP))=a+b. log(Nominal), où B sera le poids à utiliser (arrondi au nombre entier le plus proche)
b) Ajustez les données de chaque cycle à un modèle linéaire (PAR = a+b Nominal) et un modèle quadratique (PAR = a+B Nominal+c Nominal^2)
c) Calculez à rebours la concentration trouvée pour chaque étalon et comparez-la à la concentration nominale pour obtenir le biais.
d) Évaluer le biais sur toute la gamme d'étalonnage et choisir le modèle en fonction du biais

Cette question essaie de faire c). Les messages sur la liste de diffusion R suggèrent qu'il n'est pas approprié de simplement faire la régression avec les termes inversés, je peux faire manuellement le calcul pour le modèle linéaire, mais j'ai du mal avec le modèle quadratique. En cherchant sur la liste de diffusion R, il semble que d'autres personnes veulent faire la même chose.

Answer 1

OK, j'ai dû essayer ça, après avoir regardé plusieurs choses, j'ai écrit une fonction pour trouver les racines d'une équation quadratique.

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
}

donc, étant donné les données originales

> 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

nous pouvons faire l'analyse, trouver les co-efficients et ensuite trouver l'inverse, en utilisant des limites raisonnables pour les valeurs nominales que nous attendons...

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

J'espère que cela vous aidera....

Answer 2

Cette question a été soulevée sur la liste R mialing (cf. aquí ). Le post suggère un paquet que je n'arrive pas à trouver, mais ces paquets pourraient être utiles ( calibre , chimie )