I have two data sets that each have the following form:
x y sigma
13 1495.00 0.07
15 1700.91 0.09
...
basically where $x$ and $y$ are given, but also the uncertainty of $y$, denoted $\sigma$, is different for each point. The uncertainties reported are claimed to be the "68% credible region" around $y$.
I want to know whether the slope of the linear regression between $x$ and $y$ differs between the two data sets. It's easy enough for me to just do (in R)
model1 <- lm(data1$y ~ data$x)
slope1 <- coef(model1)[2]
d_slope1 <- summary(model1)$coefficients[2,2]
and the same for data2 and then check to see if they are outside of each other's uncertainty. But this ignores the information of $\sigma$. How can I incorporate $\sigma$ into my calculations?
My idea has been to add rnrom(nrow(data), 0, sigma) to y some large number of times and calculate the average standard error of the slope over all of these simulated noise additions, something like:
out <- rep(0, 1000)
for (ii in 1:1000) {
noisy_y <- data1$y + rnrom(nrow(data1), 0, sigma)
out[ii] <- summary(lm(noisy_y ~ data1$x))$coefficients[2,2]
}
better_d_slope1 <- mean(out)
Is this the correct way to do it?
Additional complication: I'm actually doing a weighted fit, where the weights come from another data set. I would think that one would normally set the weights to be $1/\sigma$, but I already have weights. (Maybe I should multiply them?)
