Comment cookdistance

.vscode/settings.json

@@ -0,0 +1,2 @@
+{
+}

clusters_clt.tex

@@ -0,0 +1,24 @@
+\setlength{\LTpost}{0mm}
+\begin{longtable}{l|rrrrrrrr}
+\toprule
+\multicolumn{1}{l}{} & \multicolumn{2}{c}{\emph{G} = 50} & \multicolumn{2}{c}{\emph{G} = 100} & \multicolumn{2}{c}{\emph{G} = 150} & \multicolumn{2}{c}{\emph{G} = 200} \\ 
+\cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} \cmidrule(lr){8-9}
+\multicolumn{1}{l}{} & \(\frac{\sqrt{G}\bar{X}_{n}}{\bar{\sigma}_{n}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\hat{\sigma}_{n}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\bar{\sigma}_{n}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\hat{\sigma}_{n}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\bar{\sigma}_{n}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\hat{\sigma}_{n}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\bar{\sigma}_{n}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\hat{\sigma}_{n}}\) \\ 
+\midrule\addlinespace[2.5pt]
+\multicolumn{9}{l}{\(\alpha=1.5, \beta=1.9\)} \\ 
+\midrule\addlinespace[2.5pt]
+10\% & $0.03$ & $0.06$ & $0.04$ & $0.03$ & $0.03$ & $0.04$ & $0.05$ & $0.06$ \\ 
+5\% & $0.02$ & $0.03$ & $0.02$ & $0.01$ & $0.00$ & $0.01$ & $0.04$ & $0.06$ \\ 
+1\% & $0.01$ & $0.00$ & $0.01$ & $0.00$ & $0.00$ & $0.01$ & $0.01$ & $0.02$ \\ 
+\midrule\addlinespace[2.5pt]
+\multicolumn{9}{l}{\(\alpha=1.5, \beta=2.1\)} \\ 
+\midrule\addlinespace[2.5pt]
+10\% & $0.05$ & $0.03$ & $0.03$ & $0.08$ & $0.03$ & $0.03$ & $0.02$ & $0.08$ \\ 
+5\% & $0.02$ & $0.02$ & $0.01$ & $0.05$ & $0.02$ & $0.01$ & $0.02$ & $0.05$ \\ 
+1\% & $0.01$ & $0.00$ & $0.00$ & $0.02$ & $0.00$ & $0.00$ & $0.00$ & $0.01$ \\ 
+\bottomrule
+\end{longtable}
+\begin{minipage}{\linewidth}
+The table shows the rejection rates for \(G=\{50,100,150,200\}\) and \(\alpha=1.5\) and \(\beta=1.9\) and \(\beta=2.1\). The Monte Carlo is based on 10,000 simulations. The simulation standard errors are: 0.009 for \(\alpha=10\%\), 0.007 for \(\alpha=5\%\), and 0.0031 for \(\alpha=1\%\).\\
+\end{minipage}
+

row2col2.tex

@@ -0,0 +1,15 @@
+\begin{tabular}{cccccccc}
+\toprule
+\multirow{2}{*}{level0} & \multirow{2}{*}{level1}&\multicolumn{3}{c}{col01}&\multicolumn{3}{c}{col02}\tabularnewline
+\cmidrule(lr){3-5}
+\cmidrule(lr){6-8}
+&&c1&c2&c3&c1&c2&c3\tabularnewline
+\midrule
+\midrule
+\multirow{2}{*}{r01}&r1&1.80e+00& (1.88e+00, -1.01e+00)& (1.27e+00, -1.31e+00)&1.80e+00& (1.88e+00, -1.01e+00)& (1.27e+00, -1.31e+00)\tabularnewline
+&r2&9.25e-01& (6.62e-01, -2.21e-01)& (-3.72e-01, 2.65e+00)&9.25e-01& (6.62e-01, -2.21e-01)& (-3.72e-01, 2.65e+00)\tabularnewline
+\midrule
+\multirow{2}{*}{r02}&r1&1.80e+00& (1.88e+00, -1.01e+00)& (1.27e+00, -1.31e+00)&1.80e+00& (1.88e+00, -1.01e+00)& (1.27e+00, -1.31e+00)\tabularnewline
+&r2&9.25e-01& (6.62e-01, -2.21e-01)& (-3.72e-01, 2.65e+00)&9.25e-01& (6.62e-01, -2.21e-01)& (-3.72e-01, 2.65e+00)\tabularnewline
+\bottomrule
+\end{tabular}

src/glmtools.jl

@@ -557,17 +557,17 @@ loglik_apweights_obs(::Poisson, y, μ, wt, ϕ, sumwt, n) = wt*logpdf(Poisson(μ)
 loglik_apweights_obs(d::NegativeBinomial, y, μ, wt, ϕ, sumwt, n) = wt*logpdf(NegativeBinomial(d.r, d.r/(μ+d.r)), y)
 
 ## Studentized Pearson residuals
-function pearson_residuals(m::GeneralizedLinearModel)
-    r = m.rr
-    wts = r.wts
-    y, η, μ = r.y, r.eta, r.mu
-    h = leverage(m)
-    sqrt(dispersion(m)).*((y .- μ).*sqrt.(wts)) ./ sqrt.(dispersion(m).*glmvar.(r.d, μ))
-end
-
-function ccooksdistance(m::GeneralizedLinearModel)
-    h = leverage(m)
-    hh = h./(1.0.-h).^2
-    Rp = pearson_residuals(m)
-    (Rp.^2).*hh./(dispersion(m)^2*dof(m))
-end
+# function pearson_residuals(m::GeneralizedLinearModel)
+#     r = m.rr
+#     wts = r.wts
+#     y, η, μ = r.y, r.eta, r.mu
+#     h = leverage(m)
+#     sqrt(dispersion(m)).*((y .- μ).*sqrt.(wts)) ./ sqrt.(dispersion(m).*glmvar.(r.d, μ))
+# end
+
+# function cooksdistance(m::GeneralizedLinearModel)
+#     h = leverage(m)
+#     hh = h./(1.0.-h).^2
+#     Rp = pearson_residuals(m)
+#     (Rp.^2).*hh./(dispersion(m)^2*dof(m))
+# end

test.tex

@@ -0,0 +1,20 @@
+\begin{longtable}{l|rrrrrr}
+\toprule
+\multicolumn{1}{l}{} & \multicolumn{2}{c}{\emph{G} = 50} & \multicolumn{2}{c}{\emph{G} = 100} & \multicolumn{2}{c}{\emph{G} = 150} \\ 
+\cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7}
+\multicolumn{1}{l}{} & \(\frac{\sqrt{G}\bar{X}_{n}}{\bar{\sigma}_{n}^{2}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\bar{\sigma}_{n}^{2}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\bar{\sigma}_{n}^{2}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\bar{\sigma}_{n}^{2}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\bar{\sigma}_{n}^{2}}\) & \(\frac{\sqrt{G}\bar{X}_{n}}{\bar{\sigma}_{n}^{2}}\) \\ 
+\midrule\addlinespace[2.5pt]
+\multicolumn{7}{l}{α=1.5, β=1.9} \\ 
+\midrule\addlinespace[2.5pt]
+10\% & $0.06$ & $0.07$ & $0.03$ & $0.03$ & $0.06$ & $0.07$ \\ 
+5\% & $0.03$ & $0.03$ & $0.01$ & $0.01$ & $0.03$ & $0.04$ \\ 
+1\% & $0.00$ & $0.00$ & $0.00$ & $0.01$ & $0.01$ & $0.01$ \\ 
+\midrule\addlinespace[2.5pt]
+\multicolumn{7}{l}{α=1.5, β=2.1} \\ 
+\midrule\addlinespace[2.5pt]
+10\% & $0.03$ & $0.05$ & $0.06$ & $0.07$ & $0.03$ & $0.04$ \\ 
+5\% & $0.01$ & $0.02$ & $0.02$ & $0.02$ & $0.01$ & $0.04$ \\ 
+1\% & $0.00$ & $0.00$ & $0.00$ & $0.01$ & $0.00$ & $0.01$ \\ 
+\bottomrule
+\end{longtable}
+

test/junk.jl

@@ -0,0 +1,30 @@
+rng = StableRNG(123)
+df = DataFrame(x_1 = randn(rng, 10), x_2 = randn(rng, 10), y = randn(rng, 10), )
+df.xx_1 = df.x_1
+df.xx_2 = df.x_2
+df.d = rand(rng, 0:1, 10)
+df.w = rand(rng, 10)
+frm0 = @formula(y ~ x_1 + x_2)
+frm1 = @formula(y ~ x_1 + xx_2 + + x_2 + xx_1)
+frmp0 = @formula(d ~ x_1 + x_2)
+frmp1 = @formula(d ~ x_1 + xx_2 + + x_2 + xx_1)
+
+probit0 = glm(frmp0, df, Binomial(), ProbitLink(), method=:qr, wts=df.w)
+
+W = Diagonal(sqrt.(probit0.rr.wrkwt))
+X = modelmatrix(probit0, weighted=false)
+
+lm0 = lm(frm0, df, wts=df.w)
+X = modelmatrix(probit0, weighted=false)
+W = Diagonal(sqrt.(probit0.rr.wts))
+v1 = diag(W*X*inv(X'W*W*X)*X'W)
+v2 = diag(X*inv(X'W*W*X)*X').*probit0.rr.wts
+
+pp = lm0.pp
+rnk = rank(pp.chol)
+
+X/pp.chol
+
+p = pp.chol.p[1:rnk]
+sum(x->x^2, view(X, :, p)/view(pp.chol.U, p, p), dims=2)
+