Last updated: 2018-05-05
workflowr checks: (Click a bullet for more information) ✖ R Markdown file: uncommitted changes
The R Markdown is untracked by Git. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.
✖ Environment: objects present
wflow_publish or wflow_build to ensure that the code is always run in an empty environment. The following objects were defined in the global environment when these results were created:
| Name | Class | Size |
|---|---|---|
| k | integer | 48 bytes |
| m | numeric | 2 Kb |
| simu.out | list | 28.7 Kb |
| simu_study | function | 14.3 Kb |
| smash.gen | function | 122.3 Kb |
✔ Seed:
set.seed(20180501)
The command set.seed(20180501) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
✔ Session information: recorded
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
✔ Repository version: b7e89a3
wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: docs/figure/
Ignored: log/
Untracked files:
Untracked: analysis/nugget.Rmd
Untracked: analysis/poireview.Rmd
Untracked: analysis/simu1.Rmd
Unstaged changes:
Modified: analysis/index.Rmd
A random variable \(X\) has a Poisson distribution with parameter \(\mu\) if it takes integer values \(y = 0, 1, 2, \dots\) with probability \(P(X=x)=\frac{e^{-\mu}\mu^x}{x!}\) where \(\mu>0\). The mean and variance of \(X\) is \(E(X)=Var(X)=\mu\).
\(f_Y(y;\theta,\phi)=\exp[(y\theta-b(\theta))/a(\phi)+c(y,\phi)]\). If \(\phi\) is known, then \(\theta\) is called canonical parameter. For example, \(\theta=\mu, \phi=\sigma^2\) in normal distribution.
\(E(Y)=\mu=b'(\theta)\), \(Var(Y)=b''(\theta)a(\phi)\). \(b''(\theta)\) is variance function and is denoted \(V(\mu)\) if it’s a function of \(\mu\). \(a(\phi)\) has the form \(a(\phi)=\frac{\phi}{a}\), where \(\phi\) is the dispersion parameter and is constant over observations, and \(a\) is the weight on each observation.
For Poisson distribution, \(\phi=1, b(\theta)=\exp(\theta),V(\mu)=1\)
\(y_i=x_i\beta+\epsilon_i\) where \(Var(\epsilon_i|x_i)=f(x_i)\neq constant\) and \(Cov(\epsilon_i,\epsilon_j)=0\).
Let \(w_i\propto 1/\sigma_i^2\), we have \(E(W^{-1/2}y|X)=W^{-1/2}X\beta\) and \(\hat\beta=(X^TWX)^{-1}X^TWy\).
Weight: \(W=V^{-1}(\frac{d\mu}{d\eta})^2\)
Dependent variate: \(z=\eta+(y-\mu)\frac{d\eta}{d\mu}\)
\(z\) is a linearized form of link function applied to the data: \(g(y)\simeq g(\mu)+(y-\mu)g'(\mu)=\eta+(y-\mu)\frac{d\eta}{d\mu}\).
Derivation:
The score function of glm model is \[s(\beta_j)=\frac{\partial l}{\partial \beta_j}=\Sigma_{i=1}^na(\phi_i)^{-1}V(\mu_i)^{-1}(y_i-\mu_i)x_{ij}\frac{d\eta_i}{d\mu_i}\] \[=\Sigma_{i=1}^na_iV(\mu_i)^{-1}(y_i-\mu_i)x_{ij}\frac{d\eta_i}{d\mu_i}\]
A method of solving score equations is the iterative algorithm Fisher’s Method of Scoring.
It turns out that the updates can be written in the form of weighted least squares where the weight matrix is \(W\) defined above and the variable \(z\) is called adjusted variable(working variable).
\(\log(\mu)=X\beta, \eta=\log(\mu)\) and \(z=\log(\mu)+\frac{y-\mu}{\mu}\)
\(z_i=\log(\hat\mu_i)+\frac{y_i-\hat\mu_i}{\hat\mu_i}\) and \(w_{ii}=\hat\mu_i\)
Notice that \(E(z_i)=\log(\mu_i)\) and \(Var(z_i)=\frac{1}{\mu_i^2}Var(y_i)=1/\mu_i\).
Gaussian nonparametric regression(Gaussian sequence model) is defined as \(y_i=\mu_i+\sigma_iz_i\) where \(z_i\sim N(0,1)\), \(i=1,2,\dots,T\). In multivariate form, it can be formulated as \(y|\mu\sim N_T(\mu,D)\) where D is the diagonal matrix with diagonal elements (\(\sigma_1^2,\dots,\sigma_T^2\)). Applying a discrete wavelet transform(DWT) represented as an orthogonal matrix \(W\), we have \(Wy|W\mu\sim N(W\mu,WDW^T)\) which is \(\tilde y|\tilde \mu\sim N(\tilde\mu,WDW^T)\). If \(\mu\) has spatial structure, then \(\tilde\mu\) would be sparse.
In heterokedastic variance case, we only use the diagonal of \(WDW^T\) so we can apply EB shrinkage to \(\tilde y_j|\tilde\mu_j\sim N(\tilde\mu_j,w_j^2)\) where \(w_j^2=\Sigma_{t=1}^T\sigma_t^2W_{jt}^2\).
If \(D\) is unknown, we estimate it using shrinkage methods under the assumption that the variances are also spatially structured.
Define \(Z_t^2=(y_t-\mu_t)^2\). Since \(y_t-\mu_t\sim N(0,\sigma_t^2)\), \(Z_t^2\sim\sigma_t^2\chi_1^2\) and \(E(Z_t^2)=\sigma_t^2\). Though \(Z_t^2\) is chi-squared distributed, we treat it as Gaussian. We use \(\frac{2}{3}Z_t^4\) to estimate \(Var(Z_t^2)\).(\(Var(Z_t^2)=2\sigma_t^4\), \(E(Z_t^4)=3\sigma_t^4\).)
Suppose \(Y_t = \mu_t + N(0,s^2_t) + N(0,\sigma^2)\) where \(s^2_t\) is known and the mean vector \(\mu\) and (constant) \(\sigma\) are to be estimated.
If assume \(\sigma\) is known too, then it is equivalent to the above heterokedastic variance case.
For Poisson data, \(X_t\sim Poi(m_t)\), let \(Y_t=\log(m_t)+\frac{x_t-m_t}{m_t}\) and set \(\mu_t=\log(m_t)\), \(s_t^2=\frac{1}{m_t}\). We can then estimate \(\mu_t\) and \(\sigma^2\) using the above generalized smash model.
Algorithm:
Some Rational: consider the Taylor series expansion of \(\log(X_t)\) around \(\lambda_t\): \(\log(X_t)=\log(\lambda_t)+\frac{1}{\lambda_t}(X_t-\lambda_t)-\frac{1}{2\lambda_t^2}(X_t-\lambda_t)^2+\dots\)
sessionInfo()R version 3.4.0 (2017-04-21)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] smashr_1.1-1 workflowr_1.0.1 rmarkdown_1.8
loaded via a namespace (and not attached):
[1] Rcpp_0.12.16 highr_0.6 compiler_3.4.0
[4] git2r_0.21.0 R.methodsS3_1.7.1 R.utils_2.6.0
[7] bitops_1.0-6 iterators_1.0.8 tools_3.4.0
[10] digest_0.6.13 evaluate_0.10 lattice_0.20-35
[13] Matrix_1.2-9 foreach_1.4.3 rstudioapi_0.7
[16] yaml_2.1.19 parallel_3.4.0 stringr_1.3.0
[19] knitr_1.20 REBayes_1.3 caTools_1.17.1
[22] rprojroot_1.3-2 grid_3.4.0 glue_1.2.0
[25] data.table_1.10.4-3 R6_2.2.2 ashr_2.2-7
[28] magrittr_1.5 whisker_0.3-2 backports_1.0.5
[31] codetools_0.2-15 htmltools_0.3.5 MASS_7.3-47
[34] assertthat_0.2.0 wavethresh_4.6.8 stringi_1.1.6
[37] Rmosek_8.0.69 doParallel_1.0.11 pscl_1.4.9
[40] truncnorm_1.0-7 SQUAREM_2017.10-1 crayon_1.3.4
[43] R.oo_1.21.0 This reproducible R Markdown analysis was created with workflowr 1.0.1