Regresi ketika setiap titik memiliki ketidakpastiannya sendiri dalam

12

Saya membuat pengukuran dua variabel dan . Mereka berdua mengetahui ketidakpastian dan terkait dengannya. Saya ingin menemukan hubungan antara dan . Bagaimana saya bisa melakukannya?x y σ x σ y x ynxyσxσyxy

EDIT : setiap memiliki terkait dengannya, dan sama dengan .σ x , i y ixiσx,iyi


Contoh R yang dapat direproduksi:

## pick some real x and y values 
true_x <- 1:100
true_y <- 2*true_x+1

## pick the uncertainty on them
sigma_x <- runif(length(true_x), 1, 10) # 10
sigma_y <- runif(length(true_y), 1, 15) # 15

## perturb both x and y with noise 
noisy_x <- rnorm(length(true_x), true_x, sigma_x)
noisy_y <- rnorm(length(true_y), true_y, sigma_y)

## make a plot 
plot(NA, xlab="x", ylab="y",
    xlim=range(noisy_x-sigma_x, noisy_x+sigma_x), 
    ylim=range(noisy_y-sigma_y, noisy_y+sigma_y))
arrows(noisy_x, noisy_y-sigma_y, 
       noisy_x, noisy_y+sigma_y, 
       length=0, angle=90, code=3, col="darkgray")
arrows(noisy_x-sigma_x, noisy_y,
       noisy_x+sigma_x, noisy_y,
       length=0, angle=90, code=3, col="darkgray")
points(noisy_y ~ noisy_x)

## fit a line 
mdl <- lm(noisy_y ~ noisy_x)
abline(mdl)

## show confidence interval around line 
newXs <- seq(-100, 200, 1)
prd <- predict(mdl, newdata=data.frame(noisy_x=newXs), 
    interval=c('confidence'), level=0.99, type='response')
lines(newXs, prd[,2], col='black', lty=3)
lines(newXs, prd[,3], col='black', lty=3)

regresi linier tanpa mempertimbangkan kesalahan dalam variabel

Masalah dengan contoh ini adalah saya pikir itu mengasumsikan bahwa tidak ada ketidakpastian dalam . Bagaimana saya bisa memperbaikinya?x

rhombidodecahedron
sumber
Benar, lmcocok dengan model regresi linier, yaitu: model ekspektasi terhadap , di mana jelas adalah acak dan dianggap dikenal. Untuk menghadapi ketidakpastian dalam Anda membutuhkan model yang berbeda. P ( Y | X ) Y X XYP(Y|X)YXX
conjugateprior
1
Untuk kasus Anda yang agak khusus (univariat dengan rasio tingkat kebisingan yang diketahui untuk X dan Y) Deming regression akan melakukan trik, misalnya Demingfungsi dalam paket R MethComp .
conjugateprior
1
@conjugateprior Terima kasih, ini terlihat menjanjikan. Saya bertanya-tanya: apakah regresi Deming masih berfungsi jika saya memiliki varian yang berbeda (tetapi masih diketahui) pada masing-masing individu x dan y? yaitu jika x adalah panjang, dan saya menggunakan penggaris dengan presisi yang berbeda untuk mendapatkan masing
rhombidodecahedron
Saya pikir mungkin cara untuk menyelesaikannya ketika ada varian yang berbeda untuk setiap pengukuran menggunakan metode York. Adakah yang tahu kalau ada implementasi R dari metode ini?
rhombidodecahedron
1
@rhombidodecahedron Lihat "dengan kesalahan yang diukur" sesuai dengan jawaban saya di sana: stats.stackexchange.com/questions/174533/… (yang diambil dari dokumentasi deming paket).
Roland

Jawaban:

9

Lθγ

(x,y):cos(θ)x+sin(θ)y=γ.

(x,y)

d(x,y;L)=cos(θ)x+sin(θ)yγ.

xiσi2yiτi2xiyi

Var(d(xi,yi;L))=cos2(θ)σi2+sin2(θ)τi2.

θγ

σiτi0


τiσixn=8

Angka

Garis yang benar ditampilkan dalam titik-titik biru. Sepanjang itu titik-titik aslinya diplot sebagai lingkaran kosong. Panah abu-abu menghubungkannya ke titik yang diamati, diplot sebagai cakram hitam pekat. Solusinya digambarkan sebagai garis merah solid. Meskipun ada penyimpangan besar antara nilai-nilai yang diamati dan aktual, solusinya sangat dekat dengan garis yang benar dalam wilayah ini.

#
# Generate data.
#
theta <- c(1, -2, 3) # The line is theta %*% c(x,y,-1) == 0
theta[-3] <- theta[-3]/sqrt(crossprod(theta[-3]))
n <- 8
set.seed(17)
sigma <- rexp(n, 1/2)
tau <- rexp(n, 1)
u <- 1:n
xy.0 <- t(outer(c(-theta[2], theta[1]), 0:(n-1)) + c(theta[3]/theta[1], 0))
xy <- xy.0 + cbind(rnorm(n, sd=sigma), rnorm(n, sd=tau))
#
# Fit a line.
#
x <- xy[, 1]
y <- xy[, 2]
f <- function(phi) { # Negative log likelihood, up to an additive constant
  a <- phi[1]
  gamma <- phi[2]
  sum((x*cos(a) + y*sin(a) - gamma)^2 / ((sigma*cos(a))^2 + (tau*sin(a))^2))/2
}
fit <- lm(y ~ x) # Yields starting estimates
slope <- coef(fit)[2]
theta.0 <- atan2(1, -slope)
gamma.0 <- coef(fit)[1] / sqrt(1 + slope^2)
sol <- nlm(f,c(theta.0, gamma.0))
#
# Plot the data and the fit.
#
theta.hat <- sol$estimate[1] %% (2*pi)
gamma.hat <- sol$estimate[2]
plot(rbind(xy.0, xy), type="n", xlab="x", ylab="y")
invisible(sapply(1:n, function(i) 
  arrows(xy.0[i,1], xy.0[i,2], xy[i,1], xy[i,2], 
         length=0.15, angle=20, col="Gray")))
points(xy.0)
points(xy, pch=16)
abline(c(theta[3] / theta[2], -theta[1]/theta[2]), col="Blue", lwd=2, lty=3)
abline(c(gamma.hat / sin(theta.hat), -1/tan(theta.hat)), col="Red", lwd=2)
whuber
sumber
+1. Sejauh yang saya mengerti, ini juga menjawab Q yang lebih lama ini: stats.stackexchange.com/questions/178727 ? Kita harus menutupnya sebagai duplikat.
Amuba kata Reinstate Monica
Juga, sesuai komentar saya untuk jawaban di utas itu, sepertinya demingfungsi dapat menangani kesalahan variabel juga. Mungkin harus menghasilkan cocok sangat mirip dengan milikmu.
Amuba kata Reinstate Monica
Saya bertanya-tanya apakah alur diskusi lebih masuk akal jika Anda mengganti tempat 2 paragraf di atas & di bawah angka?
gung - Reinstate Monica
3
Saya diingatkan pagi ini (oleh seorang pemilih) bahwa pertanyaan ini telah ditanyakan dan dijawab dengan berbagai cara, dengan kode kerja, beberapa tahun yang lalu di situs Mathematica SE .
Whuber
Apakah solusi ini punya nama? dan mungkin sumber daya untuk bacaan lebih lanjut (selain situs SE Mathematica yang saya maksud)?
JustGettin Mulai
0

Optimalisasi kemungkinan maksimum untuk kasus ketidakpastian x dan y telah ditangani oleh York (2004). Berikut adalah kode R untuk fungsinya.

"YorkFit", ditulis oleh Rick Wehr, 2011, diterjemahkan ke dalam R oleh Rachel Chang

Rutin universal untuk menemukan garis lurus terbaik sesuai dengan data dengan variabel, kesalahan berkorelasi, termasuk estimasi kesalahan dan goodness of fit, mengikuti Persamaan. (13) dari York 2004, American Journal of Physics, yang didasarkan pada York 1969, Earth and Planetary Sciences Letters

YorkFit <- fungsi (X, Y, Xstd, Ystd, Ri = 0, b0 = 0, printCoefs = 0, makeLine = 0, eps = 1e-7)

X, Y, Xstd, Ystd: gelombang yang berisi titik X, titik Y, dan deviasi standar mereka

PERINGATAN: Xstd dan Ystd tidak boleh nol karena ini akan menyebabkan Xw atau Yw menjadi NaN. Gunakan nilai yang sangat kecil sebagai gantinya.

Ri: koefisien korelasi untuk kesalahan X dan Y - panjang 1 atau panjang X dan Y

b0: tebakan awal kasar untuk kemiringan (dapat diperoleh dari standar kuadrat-terkecil tanpa kesalahan)

printCoefs: set sama dengan 1 untuk menampilkan hasil di jendela perintah

makeLine: set sama dengan 1 untuk menghasilkan gelombang Y untuk garis fit

Mengembalikan matriks dengan intersep dan kemiringan ditambah ketidakpastiannya

Jika tidak ada tebakan awal untuk b0, maka gunakan saja OLS jika (b0 == 0) {b0 = lm (Y ~ X) $ koefisien [2]}

tol = abs(b0)*eps #the fit will stop iterating when the slope converges to within this value

a, b: intersep dan kemiringan akhir a.err, b.err: estimasi ketidakpastian dalam intersep dan kemiringan

# WAVE DEFINITIONS #

Xw = 1/(Xstd^2) #X weights
Yw = 1/(Ystd^2) #Y weights


# ITERATIVE CALCULATION OF SLOPE AND INTERCEPT #

b = b0
b.diff = tol + 1
while(b.diff>tol)
{
    b.old = b
    alpha.i = sqrt(Xw*Yw)
    Wi = (Xw*Yw)/((b^2)*Yw + Xw - 2*b*Ri*alpha.i)
    WiX = Wi*X
    WiY = Wi*Y
    sumWiX = sum(WiX, na.rm = TRUE)
    sumWiY = sum(WiY, na.rm = TRUE)
    sumWi = sum(Wi, na.rm = TRUE)
    Xbar = sumWiX/sumWi
    Ybar = sumWiY/sumWi
    Ui = X - Xbar
    Vi = Y - Ybar

    Bi = Wi*((Ui/Yw) + (b*Vi/Xw) - (b*Ui+Vi)*Ri/alpha.i)
    wTOPint = Bi*Wi*Vi
    wBOTint = Bi*Wi*Ui
    sumTOP = sum(wTOPint, na.rm=TRUE)
    sumBOT = sum(wBOTint, na.rm=TRUE)
    b = sumTOP/sumBOT

    b.diff = abs(b-b.old)
  }     

   a = Ybar - b*Xbar
   wYorkFitCoefs = c(a,b)

# ERROR CALCULATION #

Xadj = Xbar + Bi
WiXadj = Wi*Xadj
sumWiXadj = sum(WiXadj, na.rm=TRUE)
Xadjbar = sumWiXadj/sumWi
Uadj = Xadj - Xadjbar
wErrorTerm = Wi*Uadj*Uadj
errorSum = sum(wErrorTerm, na.rm=TRUE)
b.err = sqrt(1/errorSum)
a.err = sqrt((1/sumWi) + (Xadjbar^2)*(b.err^2))
wYorkFitErrors = c(a.err,b.err)

# GOODNESS OF FIT CALCULATION #
lgth = length(X)
wSint = Wi*(Y - b*X - a)^2
sumSint = sum(wSint, na.rm=TRUE)
wYorkGOF = c(sumSint/(lgth-2),sqrt(2/(lgth-2))) #GOF (should equal 1 if assumptions are valid), #standard error in GOF

# OPTIONAL OUTPUTS #

if(printCoefs==1)
 {
    print(paste("intercept = ", a, " +/- ", a.err, sep=""))
    print(paste("slope = ", b, " +/- ", b.err, sep=""))
  }
if(makeLine==1)
 {
    wYorkFitLine = a + b*X
  }
 ans=rbind(c(a,a.err),c(b, b.err)); dimnames(ans)=list(c("Int","Slope"),c("Value","Sigma"))
return(ans)
 }
Steven Wofsy
sumber
Juga perhatikan, paket R "IsoplotR" termasuk fungsi york (), memberikan hasil yang sama dengan kode YorkFit di sini.
Steven Wofsy