Utilizzo di incertezza sperimentale durante l'esecuzione della regressione lineare in R

Question

Ho un set di dati sperimentali che mi piacerebbe inserire in un polinomio. I dati comprendono una variabile indipendente, la variabile dipendente ed un'incertezza nella misurazione di quest'ultimo, ad esempio

2000 0.2084272 0.002067834 
2500 0.207078125 0.001037248 
3000 0.2054202 0.001959138 
3500 0.203488075 0.000328942 
4000 0.2013152 0.000646088 
4500 0.198933825 0.001375657 
5000 0.196375 0.000908696 
5500 0.193668575 0.00014721 
6000 0.1908432 0.000526976 
6500 0.187926325 0.001217318 
7000 0.1849442 0.000556495 
7500 0.181921875 0.000401883 
8000 0.1788832 0.001446992 
8500 0.175850825 0.0
9000 0.1728462 0.001676249 
9500 0.169889575 0.001011735 
10000 0.167  0.000326678 

(colonne x, y, + -y).

posso effettuare un polinomio usando il suddetto con ad esempio

mydata = read.table("example.txt") 
model <- lm(V2~V1+I(V1^2)+I(V1^3)+I(V1^4), data = mydata) 

ma questo non utilizzare i valori di incertezza. Come informo R che la terza colonna dell'insieme di dati è un'incertezza e che pertanto dovrebbe essere utilizzata nell'analisi di regressione?

Answer 1

Con errore di misurazione nella variabile dipendente non correlato con le variabili indipendenti, i coefficienti stimati sono imparziali ma gli errori standard sono troppo piccoli. Qui è il riferimento che ho usato (Pagina 1 & 2): http://people.stfx.ca/tleo/econ370term2lec4.pdf

Penso che basta per regolare gli errori standard da quelli calcolati da LM(). Ed è quello che ho cercato di fare nel codice qui sotto. Non sono una persona stat, quindi potresti voler postare per convalidare e chiedere intuizioni migliori.

Per l'esempio di seguito, ho assunto che la colonna "incertezza" sia la deviazione standard (o l'errore standard). Per semplicità, ho modificato il modello solo: y ~ x.

# initialize dataset 
df <- data.frame(
     x = c(2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000,7500,8000,8500,9000,9500,10000), 
     y = c(0.2084272,0.207078125,0.2054202,0.203488075,0.2013152,0.198933825, 0.196375,0.193668575, 0.1908432, 0.187926325, 0.1849442, 0.181921875, 0.1788832, 0.175850825, 0.1728462,0.169889575, 0.167), 
     y.err = c(0.002067834, 0.001037248, 0.001959138, 0.000328942, 0.000646088, 0.001375657, 0.000908696, 0.00014721, 0.000526976, 0.001217318, 0.000556495, 0.000401883, 0.001446992, 0.0, 0.001676249, 0.001011735, 0.000326678) 
) 

df 

# model regression 
model <- lm(y~x, data = df) 
summary(model) 

# get the variance of the measurement error 
# thanks to: http://schools-wikipedia.org/wp/v/Variance.htm 
# law of total variance 
pooled.var.y.err <- mean((df$y.err)^2) + var((df$y.err)^2) 

# get variance of beta from model 
# thanks to: http://stats.stackexchange.com/questions/44838/how-are-the-standard-errors-of-coefficients-calculated-in-a-regression 
X <- cbind(1, df$x) 
#  (if you add more variables to the model you need to modify the following line) 
var.betaHat <- anova(model)[[3]][2] * solve(t(X) %*% X) 

# add betaHat variance to measurement error variance 
var.betaHat.adj <- var.betaHat + pooled.var.y.err 

# calculate adjusted standard errors 
sqrt(diag(var.betaHat.adj)) 

# compare to un-adjusted standard errors 
sqrt(diag(var.betaHat))