[R-meta] FW: Welcome to the "R-sig-meta-analysis" mailing list (Digest mode) _ question

Dear Tetyana

Have you looked at the Riley et al paper referenced in the 
documentation? I think it has the formula although I have not time to 
check that now.

Michael

Dear All,

I am wondering if you can help me understanding the underlaying
calculations/formula for the confidence intervals (both) using the
predic.rma function.

I will try to explain the situation using an example below.

test_TST = rma(yi=TST_Mean, sei=TST_E, method="DL",

mods=~Age_Mean+Sex__1+Night_coded, data=moderators[-c(5,69,117,168), ])

test_TST

Mixed-Effects Model (k = 128; tau^2 estimator: DL)

tau^2 (estimated amount of residual heterogeneity):  439.3101 (SE =
263.3751) tau (square root of estimated tau^2 value):             20.9597
I^2 (residual heterogeneity / unaccounted variability): 91.95%
H^2 (unaccounted variability / sampling variability):   12.42
R^2 (amount of heterogeneity accounted for):            70.67%

Test for Residual Heterogeneity:
QE(df = 124) = 1539.9995, p-val < .0001

Test of Moderators (coefficient(s) 2,3,4):
QM(df = 3) = 183.5661, p-val < .0001

Model Results:

  		estimate       se     	zval    pval     ci.lb
ci.ub
  intrcpt      	375.7589  10.4376  36.0005  <.0001  355.3016  396.2162  ***
Age_Mean      	-1.0134     0.1345    -7.5366  <.0001   -1.2770   -0.7499
***
Sex__1           	0.0291       0.0667    0.4365  0.6624   -0.1017
0.1599
Night_coded     38.3012     4.5223     8.4694  <.0001   29.4376   47.1647
***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

I am trying to predict the estimates for the outcome, TST, based on certain
parameters using the predic.rma function:

predict.rma (test_TST, newmods=cbind(c(65),c(100),c(1)))

      pred     se    		ci.lb    ci.ub    cr.lb    	cr.ub
351.1013 5.3417 340.6318 361.5708 308.7079 393.4947


My question is how I can use estimates from the rm () model results above
(the formula) to get the ci.lb    ci.ub    cr.lb    cr.ub manually?

To get pred value, we used the following formula which makes sense:
Intercept_estimate + Age_coeff_estimate X Age (65) +Sex_coeff_estimate X 0
[if women] or 100 [if men] +  Sleep_Night_coeff_estimate X 0 [first night]
or 1 [second night].

Which formula we should use to get in the same way:  ci.lb    ci.ub    cr.lb
cr.ub?

Thank you in advance and sorry if it sounds somewhat confusing,
Tetyana

_______________________________________________
R-sig-meta-analysis mailing list
R-sig-meta-analysis at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis

---
This email has been checked for viruses by AVG.
https://www.avg.com

Michael
http://www.dewey.myzen.co.uk/home.html

Hi Tetyana,

Let V_b be the estimated variance-covariance matrix of the fixed effects, which is computed with:

V_b = (X'WX)^(-1),

where X is the model matrix (with a column of 1s for the intercept and columns for the other moderators in the model) and W = M^(-1) is the inverse of the model-implied estimated variance-covariance matrix of the outcomes (M). For this model, M is just equal to a diagonal matrix with TST_E^2 + hat(tau)^2 along the diagonal, where hat(tau)^2 is the estimated value of tau^2. And hence W is diagonal with 1/ (TST_E^2 + hat(tau)^2) along the diagonal.

You can extract X with model.matrix(test_TST), M with vcov(test_TST, type="obs"), and W with weights(test_TST, type="matrix"). See help(vcov.rma) and help(weights.rma.uni). 

Now let x be a row vector with the desired values for the moderators (again, with a 1 in the first column). Then the predicted value is given by xb, where b is the column vector with the estimated model coefficients (which you can get with test_TRT$b or cbind(coef(test_LRT))). And then:

Var[xb] = x V_b x'

and SE[xb] = sqrt(Var[xb]), so that is how you will get the SE that is given by predict() and then the lower and upper 95% CI bounds are given by

xb +- 1.96 * SE[xb]

(or rather qnorm(.975) instead of 1.96 to be precise). The bounds of the credibility/prediction interval are given by

xb +- 1.96 * sqrt(Var[xb] + hat(tau)^2).

P.S.: newmods=cbind(c(65),c(100),c(1)) is unnecessarily complex. This will do: newmods=c(65,100,1).

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On Behalf Of Tetyana Kendzerska
Sent: Wednesday, 18 September, 2019 5:41
To: r-sig-meta-analysis at r-project.org
Subject: [R-meta] FW: Welcome to the "R-sig-meta-analysis" mailing list (Digest mode) _ question

Dear All,

I am wondering if you can help me understanding the underlaying
calculations/formula for the confidence intervals (both) using the
predic.rma function.

I will try to explain the situation using an example below.

test_TST = rma(yi=TST_Mean, sei=TST_E, method="DL",

mods=~Age_Mean+Sex__1+Night_coded, data=moderators[-c(5,69,117,168), ])

test_TST

Mixed-Effects Model (k = 128; tau^2 estimator: DL)

tau^2 (estimated amount of residual heterogeneity):  439.3101 (SE =
263.3751) tau (square root of estimated tau^2 value):             20.9597
I^2 (residual heterogeneity / unaccounted variability): 91.95%
H^2 (unaccounted variability / sampling variability):   12.42
R^2 (amount of heterogeneity accounted for):            70.67%

Test for Residual Heterogeneity: 
QE(df = 124) = 1539.9995, p-val < .0001

Test of Moderators (coefficient(s) 2,3,4): 
QM(df = 3) = 183.5661, p-val < .0001

Model Results:

 		estimate       se     	zval    pval     ci.lb
ci.ub 
 intrcpt      	375.7589  10.4376  36.0005  <.0001  355.3016  396.2162  *** 
Age_Mean      	-1.0134     0.1345    -7.5366  <.0001   -1.2770   -0.7499
***
Sex__1           	0.0291       0.0667    0.4365  0.6624   -0.1017
0.1599 
Night_coded     38.3012     4.5223     8.4694  <.0001   29.4376   47.1647
***

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

I am trying to predict the estimates for the outcome, TST, based on certain
parameters using the predic.rma function:

predict.rma (test_TST, newmods=cbind(c(65),c(100),c(1)))

pred     se    		ci.lb    ci.ub    cr.lb    	cr.ub
351.1013 5.3417 340.6318 361.5708 308.7079 393.4947 

My question is how I can use estimates from the rm () model results above
(the formula) to get the ci.lb    ci.ub    cr.lb    cr.ub manually?  

To get pred value, we used the following formula which makes sense: 
Intercept_estimate + Age_coeff_estimate X Age (65) +Sex_coeff_estimate X 0
[if women] or 100 [if men] +  Sleep_Night_coeff_estimate X 0 [first night]
or 1 [second night].

Which formula we should use to get in the same way:  ci.lb    ci.ub    cr.lb
cr.ub?

Thank you in advance and sorry if it sounds somewhat confusing,
Tetyana