I fit your data with AD Model Builder. The R code has the right std devs (almost) but poor parameter estimates. The C code has good parameter estimates but the wrong std devs You need to multiply by sqrt(9/7) to convert from fisher information to nls stdevs
index name value std.dev
1 N 4.7079e+00 1.7120e-01
2 rate 4.7979e-01 1.9154e-02
Assume regression model:
$$ y_i = \mathbf{x}'_i \boldsymbol{\beta} + \epsilon_i$$
Let $\boldsymbol{\epsilon}$ be your vector of error terms.
If you know that $ \mathrm{Var}\left(\boldsymbol{\epsilon} \right) = \Omega$, for example:
$$ \Omega = \begin{bmatrix} \sigma^2_1 & 0 & 0& \ldots&0\\0 &\sigma^2_2 & 0 & \ldots&0\\ 0 & 0 &\sigma^2_3& \ldots & 0\\\ldots&\ldots&\ldots &\ldots&0\\0&0&0&0&\sigma^2_n\end{bmatrix} $$
Then you can more efficiently estimate $\boldsymbol{\beta}$ using generalized least squares.
The estimator $\hat{\mathbf{b}}$ for GLS is given by:
$$\hat{\mathbf{b}} = \left(X'\Omega^{-1} X \right)^{-1}\left(X'\Omega^{-1} \mathbf{y} \right) $$
The basic idea with GLS is to give observations that are more precisely observed higher weight. The danger of this approach of course is that if $\Omega$ is not correct, you can end up with something far worse than the equal weighting of regular OLS.
Note also that weighted least squares (as what would occur for this $\Omega$), is a special case of GLS.
If you just want to use OLS estimation but calculate standard errors assuming you know $\mathrm{Var}\left(\boldsymbol{\epsilon}\right)$
\begin{align*}
\mathrm{Var}\left( \hat{\mathbf{b}}_{OLS} \right) &= \mathrm{Var}\left( \left(X'X \right)^{-1}X'\left(X\mathbf{\beta} + \boldsymbol{\epsilon} \right) \right)\\
&=\left(X'X \right)^{-1}X' \mathrm{Var}\left( \boldsymbol{\epsilon} \right) X \left(X'X \right)^{-1} \\
&=\left(X'X \right)^{-1}X' \Omega X \left(X'X \right)^{-1}
\end{align*}
The variance of the GLS estimator is given by:
\begin{align*}
\mathrm{Var}\left( \hat{\mathbf{b}}_{GLS} \right) &= \mathrm{Var}\left( \left(X'\Omega^{-1} X \right)^{-1}X'\Omega^{-1}\left(X\mathbf{\beta} + \boldsymbol{\epsilon} \right) \right)\\
&=\left(X'\Omega^{-1}X \right)^{-1}X'\Omega^{-1} \mathrm{Var}\left( \boldsymbol{\epsilon} \right) \Omega^{-1} X \left(X'\Omega^{-1}X \right)^{-1} \\
&=\left(X'\Omega^{-1}X \right)^{-1}X'\Omega^{-1} \Omega\Omega^{-1} X \left(X'\Omega^{-1}X \right)^{-1} \\
&=\left(X'\Omega^{-1}X \right)^{-1}X'\Omega^{-1} X \left(X'\Omega^{-1}X \right)^{-1}\\
&=\left(X'\Omega^{-1}X \right)^{-1}
\end{align*}
Best Answer
As I have less than 50 reputations, I'll post this as an answer. A standard way to deal with these kind of problems would be to fit a multilevel model. If $y_{ij}$ is your measurement of your outcome for unit $i$ at the $j$th wave of data collection, the model assumes that
$$ y_{ij} = \alpha_j + \beta_j x_{ij} + \epsilon_{ij}$$
(note the subscripts on the coefficients) where
$$(\alpha_j, \beta_j) \sim \mathcal N(\mu, \Sigma)$$
$$\epsilon_{ij} \sim \mathcal N(0, \Omega),$$ where $\epsilon_{ij}$ and $(\alpha_j,\beta_j)$ are independently distributed. Although the distributions do not have to be Normal, the Normal distribution is the natural choice in many applications. Also it is often assumed that $\Omega = \sigma I$, although this assumption can be easily relaxed. If you have more than one predictor in your regression, $\beta_j$ would be a vector. Note that $\Sigma$ (the covariance matrix of the regression coefficients) will contain the variance of the intercept, slope, and their covariances. You might test whether they are significantly different from zero with likelihood-ratio tests. However, you have to be careful as the null-hypothesis would lie at the boundary of the parameter space.