Analytical solution of a simple regression with fixed intercept

Question

I would like to know how to find out the analytical solution of a simple linear regression with fixed intercept = 0:

$$ s = e^{-ht}$$ $$ y = -ln(s) = h\cdot t$$

Here ist the background: I have three survival probabilities $s$ at 30, 90 and 180 days. Obviously, I have at day = 0 100% survival, so I include this observation. I know that this is contested (here) but I think in this special case it makes sense. The data I use for fitting the linear regression:

>     obs
    t    s          y
1   0 1.00 0.00000000
2  30 0.98 0.02020271
3  90 0.90 0.10536052
4 180 0.80 0.22314355

If I fit with simple regression I get this:

>     (fit1 <- lm(y~t, data=obs))

Call:
lm(formula = y ~ t, data = obs)

Coefficients:
(Intercept)            t  
  -0.008464     0.001275

This can be obtained analytically if the following function is derived:

$$f(h) = \sum (y_i - ht_i)^2$$

which gives:

$$ \frac{\sum (y_i-\bar{y})\cdot (t_i-\bar{t})}{\sum (t_i-\bar{t})^2}$$

---

UPDATE 1: This is the result of the minimization of $$f(h) = \sum (y_i - c - ht_i)^2$$. The correct result (see answers): $$ \frac{\sum (y_i\cdot t_i)}{\sum t_i^2}$$

---

The analytical results is:

yc <- with(obs,y-mean(y))
tc <- with(obs, t -  mean(t))
sum(yc*tc)/sum(tc^2)
[1] 0.001275204

The same as coefficient in the fit1. Now, if I fix intercept to intercept=0 I get this:

>     (fit2 <- lm(y~0+t, data=obs))

Call:
lm(formula = y ~ 0 + t, data = obs)

Coefficients:
   t  
0.001214  

I'm wondering how I can get an analytical solution for this. How I have to consider the fix intercept in the function $f(h)$ above?

Any idea is appreciated.

---

Here is the way I constructed the data:

set.seed(123)
# Hazard ratio
h <- 0.7
# Number of observation
n <- 50
# Model: exponential
t <- rexp(n,h)
# scale to days
t <- t*365.25
hist(t)
t <- sort(t)
# Put data into a dataframe
df0 <- data.frame(t=t)
head(df0)
# Compute probablities
df0$s <- 1 - c(1:n)/n
head(df0)
# Extract survival probabilities at 30,90 and 180 days
df0$t2 <- ceiling(df0$t/30)*30
# Select survival probablity 30, 90, 180 days
library(sqldf)
obs <- sqldf("SELECT t2 t, MAX(s) s FROM df0 WHERE t2 IN (30,90,180) GROUP BY t2")
# Add survival probability=1 at day 0
obs <- rbind(data.frame(t = 0, s = 1), obs)
# s = e^(-ht)  => y = -ln(s) = h*t
obs$y <- -log(obs$s)
plot(y~t, data=obs)
fit1 <- lm(y~t, data=obs)
abline(fit1,lty=2)
fit2 <- lm(y~0+t, data=obs)
abline(fit2,lty=2, col="red")
legend("topleft", legend=c("fit1","fit2"), col=c(1,2), lty=c(2,2))

Answer 1

This question touches on a number of existing posts but I didn't find one that related to all of it.

1. You should not include the constraint as if it were an ordinary data point, with the same uncertainty (but see point 4.)

2. ordinary regression through the origin is entirely straightforward.

   Consider $y_i = \beta x_i + \epsilon_i$

   $S = \sum_i (y_i - \beta x_i)^2$

   $dS/d\beta = \sum_i -2x_i(y_i - \beta x_i)$

   Set $=0$ and solve for $\hat{\beta}$ yields $\hat{\beta} = \sum_i x_iy_i/\sum_ix_i^2$.

   (Your notation has $t$ where I have $x$ and $h$ where I have $\beta$. Your assertion under "This can be obtained analytically if the following function is derived:" is false -- the solution you offer doesn't follow from the formula for $f(h)$)

   For additional discussion about potential issues with regression through the origin, see some of our numerous threads on the topic, including

   a. Regression without intercept: deriving $\hat{β}_1$ in least squares (no matrices)

   b. What are the uses and pitfalls of regression through the origin?

   c. Regression through the origin

   d. When forcing intercept of 0 in linear regression is acceptable/advisable

3. It's not clear that it makes sense to take logs and use ordinary regression. It depends on how the noise enters the relationship between the variables. Typically when dealing with regression that must pass through the origin the noise gets smaller as $x_i \to 0$, and so it doesn't make sense to use a model with constant noise -- you would put too much weight where the data are less certain (less informative) and too little where they're more precise.

   (Again, a number of posts mentioning the issue of how the noise term enters the equation relating $x$ to observed $y$ can be found)

4. You can use a pseudo-observation approach to impose the constraint, but since you know the line passes through (0,0) you have to make that data point have zero variance (which is the same as using weighted regression but making its weight $\to \infty$ (in practice - outside of doing it algebraically, you're limited to making it very large). This corresponds to enforcing the constraint (and yields the same parameter estimate as simply solving the equating as I did in 2. above).

Answer 2

For the regression

$ y = X \beta + \epsilon $

where $y$, $beta$ and $\epsilon$ are column vectors and $X$ is a matrix. $y$ and $X$ are known. The $\beta$, denoted $\hat{\beta}$, which minimizes the sum of the squares, $\epsilon'\epsilon$, is (the proof is standard -- google for deriving regression coefficients using matrix math):

$ \hat{\beta} = (X'X)^{-1}X'y $

In the case given in the question $\hat{\beta}$ is a scalar and all other quantities are column vectors so writing $X$ as $x$ to emphasize that it is a column vector the above reduces to:

$ \hat{\beta} = x'y / x'x $

which in R can be written

crossprod(x, y) / crossprod(x)
##             [,1]
## [1,] 0.001213874

or

sum(x*y) / sum(x*x)
## [1] 0.001213874

and gives the same answer (except possibly for numeric precision and the particular R data structure used), as:

coef(lm(y ~ x + 0))
##           x 
## 0.001213874 

Note: We used the following where x equals obs$x in the question and y equals obs$t in the question:

x <- c(0, 30, 90, 180)
y <- c(0, 0.0202027073175195, 0.105360515657826, 0.22314355131421)