Regressão quando cada ponto tem a sua própria incerteza em ambos

12

Fiz n medições de duas variáveis x e y . Ambos têm incertezas conhecidas σx e σy associadas a elas. Eu quero encontrar a relação entre x e y . Como eu posso fazer isso?

EDIT : cada xEu tem um diferente , iσx,Eu associado a ele e o mesmo com yEu .


Exemplo reproduzível de R:

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

regressão linear sem considerar erros nas variáveis

O problema com este exemplo é que acho que assume que não há incertezas em x . Como posso consertar isso?

rhombidodecahedron
fonte
É verdade que lmse encaixa em um modelo de regressão linear, ou seja: um modelo da expectativa de em relação a P ( Y | X ) , no qual Y é claramente aleatório e X é considerado conhecido. Para lidar com a incerteza em X, você precisará de um modelo diferente. YP(Y|X)YXX
conjugateprior
1
Para o seu caso bastante especial (univariado com uma taxa conhecida de níveis de ruído para X e Y), a regressão de Deming fará o truque, por exemplo, a Demingfunção no pacote R MethComp .
conjugateprior
1
@conjugateprior Obrigado, isso parece promissor. Estou me perguntando: a regressão de Deming ainda funciona se eu tiver uma variação diferente (mas ainda conhecida) em cada indivíduo xey? ou seja, se os x têm comprimentos, e eu usei réguas com diferentes precisões para obter cada x
rhombidodecahedron 15/16
Eu acho que talvez a maneira de resolvê-lo quando existem diferentes variações para cada medida seja usando o método de York. Alguém sabe se existe uma implementação R deste método?
rhombidodecahedron
1
@rhombidodecahedron Veja o ajuste "com erros medidos" na minha resposta: stats.stackexchange.com/questions/174533/… (que foi retirado da documentação de deming de pacotes).
187 Roland

Respostas:

9

Seja a reta verdadeira , dada por um ângulo θ e um valor γLθγ , o conjunto

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

A distância assinada entre qualquer ponto e esta linha é(x,y)

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

xiσi2yiτi2xiyi implica que a variação dessa distância é

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

θγ

σiτi0 : isso deve eliminar quase toda a correlação entre as estimativas dos dois parâmetros.


O método funciona tão bem com o exemplo da pergunta que a linha ajustada é quase distinguível da linha verdadeira na plotagem: elas estão dentro de uma unidade ou mais uma da outra em todos os lugares. Em vez disso, neste exemplo, oτiσixn=8

Figura

A linha verdadeira é mostrada em azul pontilhado. Ao longo, os pontos originais são plotados como círculos ocos. Setas cinza as conectam aos pontos observados, plotados como discos pretos sólidos. A solução é desenhada como uma linha vermelha sólida. Apesar da presença de grandes desvios entre os valores observados e os reais, a solução está notavelmente próxima da linha correta nessa região.

#
# 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
fonte
+1. Tanto quanto eu entendo, isso também responde a este Q mais antigo: stats.stackexchange.com/questions/178727 ? Devemos fechá-lo como duplicado então.
Ameba diz Reinstate Monica
Além disso, de acordo com o meu comentário para a resposta nesse segmento, parece que deming função também pode lidar com erros de variáveis. Provavelmente deve produzir um ajuste muito semelhante ao seu.
Ameba diz Reinstate Monica
Gostaria de saber se o fluxo da discussão faz mais sentido se você mudar os locais dos 2 parágrafos acima e abaixo da figura?
gung - Reintegrar Monica
3
Fui lembrado nesta manhã (por um eleitor) que essa pergunta havia sido feita e respondida de várias maneiras, com código de trabalho, vários anos atrás, no site do Mathematica SE .
whuber
Esta solução tem um nome? e possivelmente um recurso para leitura adicional (além do site Mathematica SE, quero dizer)?
usar o seguinte código
0

A otimização de probabilidade máxima para o caso de incertezas em xey foi abordada por York (2004). Aqui está o código R para sua função.

"YorkFit", escrito por Rick Wehr, 2011, traduzido para R por Rachel Chang

Rotina universal para encontrar o melhor ajuste de linha reta para dados com erros variáveis ​​e correlatos, incluindo estimativas de erro e qualidade do ajuste, seguindo a Eq. (13) de York 2004, American Journal of Physics, que foi baseado em York 1969, Earth and Planetary Sciences Letters

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

X, Y, Xstd, Ystd: ondas contendo pontos X, pontos Y e seus desvios padrão

AVISO: Xstd e Ystd não podem ser zero, pois isso fará com que Xw ou Yw seja NaN. Use um valor muito pequeno.

Ri: coeficientes de correlação para erros X e Y - comprimento 1 ou comprimento de X e Y

b0: estimativa inicial aproximada da inclinação (pode ser obtida de um ajuste de mínimos quadrados padrão sem erros)

printCoefs: defina igual a 1 para exibir os resultados na janela de comando

makeLine: defina igual a 1 para gerar uma onda Y para a linha de ajuste

Retorna uma matriz com a interceptação e a inclinação mais suas incertezas

Se nenhuma estimativa inicial para b0 for fornecida, use OLS se (b0 == 0) {b0 = lm (Y ~ X) $ coeficientes [2]}

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

a, b: interceptação final e inclinação a.err, b.err: incertezas estimadas em interceptação e inclinação

# 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
fonte
Observe também que o pacote R "IsoplotR" inclui a função york (), fornecendo os mesmos resultados que o código YorkFit aqui.
Steven Wofsy 23/09/19