Caculate error Odd ration in R(about baseyglm )

Question

[mmh7272]

Hello i have a problem analysis Generalized linear model.

i use glm() and bayesglm() from package in arm .

glm () function is very good working and caculate odd ration well using below function

exp(cbind(Odds_Ratio_RedvBlue=coef(fit), confint(result)))

Also,

bayesglm() function is working well .

but,! caculate error odd ration using that function

Error in profile.glm(object, which = parm, alpha = (1 - level)/4, trace = trace) :

profiling has found a better solution, so original fit had not converged

what is problem ? Plz help me

Answer 1

[Kevin Blighe]

One can quite easily obtain odds ratios and their associated confidence intervals from a Bayesian GLM. Your code, exp(cbind(Odds_Ratio_RedvBlue=coef(fit), confint(result))), however, is incorrect. The confint() function accepts a model; so, you should have fit there.

Here is a completely reproducible example:

create random data

I create 2 random data matrices and then rbind them - this helps to make both of these different, as they will be generated using different seeds. One will represent cases; the other, controls.

mat <- rbind(
  matrix(rexp(10000, rate=.1), ncol=20),
  matrix(rexp(10000, rate=.1), ncol=20))

rownames(mat) <- paste0('sample', 1:nrow(mat))
colnames(mat) <- paste0('gene', 1:ncol(mat))

add a casecontrol binary categorical variable

casecontrol <- factor(c(rep('case', nrow(mat)/2), rep('control', nrow(mat)/2)),
  levels = c('control', 'case'))
mat <- data.frame(casecontrol, mat)
mat[1:5,1:5]
        casecontrol     gene1     gene2        gene3      gene4
sample1        case 12.587667  5.273483  4.622787060 16.4662629
sample2        case  4.299452 16.972602 24.363363216  1.1530915
sample3        case  1.445478  4.427493  0.911904015 15.9480892
sample4        case  2.795976 15.170870  0.005618626  0.4279287
sample5        case  9.031451 14.965446  3.215852608  0.1271495

create a formula that includes all genes

f <- as.formula(paste0('casecontrol ~ ',
  paste(colnames(mat)[2:ncol(mat)], collapse = ' + ')))
f
casecontrol ~ gene1 + gene2 + gene3 + gene4 + gene5 + gene6 + 
    gene7 + gene8 + gene9 + gene10 + gene11 + gene12 + gene13 + 
    gene14 + gene15 + gene16 + gene17 + gene18 + gene19 + gene20

fit the formula to the data

require(arm)
model <- bayesglm(f, data = mat, family = binomial(link = 'logit'))

calculate ORs with default confidence intervals

exp(cbind(OR = coef(model), confint(model)))
Waiting for profiling to be done...
                   OR     2.5 %    97.5 %
(Intercept) 1.7419832 0.9776763 3.1350934
gene1       0.9890273 0.9759172 1.0020303
gene2       0.9940674 0.9811322 1.0070403
gene3       1.0076251 0.9948488 1.0208389
gene4       0.9863225 0.9735852 0.9987108
gene5       0.9980863 0.9851915 1.0110568
gene6       0.9912441 0.9786471 1.0036911
gene7       0.9993135 0.9868953 1.0118432
gene8       1.0071035 0.9944845 1.0201311
gene9       0.9998519 0.9872245 1.0126304
gene10      0.9984163 0.9847696 1.0121882
gene11      0.9896538 0.9763643 1.0028173
gene12      0.9956987 0.9831751 1.0082592
gene13      1.0048972 0.9935748 1.0165546
gene14      0.9923246 0.9793540 1.0052122
gene15      1.0112406 0.9979852 1.0249532
gene16      0.9895135 0.9765408 1.0023046
gene17      0.9897853 0.9763983 1.0030825
gene18      0.9959242 0.9832540 1.0086531
gene19      1.0029393 0.9904289 1.0156923
gene20      1.0011999 0.9893800 1.0132243

[mmh7272]

Oh!Kevin thank you so much ! i have a last question can add to p-value? some paper add to odd ratio table including p value so i wanna try

[Kevin Blighe]

Sure, you just need to cbind the odds ratio output to this:

summary(model)$coefficients
                 Estimate  Std. Error     z value  Pr(>|z|)
(Intercept) -0.0165666475 0.309176505 -0.05358314 0.9572673
gene1       -0.0049994409 0.006340742 -0.78846302 0.4304259
gene2        0.0097938335 0.006464123  1.51510620 0.1297455
gene3       -0.0084898872 0.006224293 -1.36399228 0.1725700
gene4       -0.0048214452 0.006874556 -0.70134647 0.4830868
gene5       -0.0077344704 0.006605418 -1.17092817 0.2416277
gene6       -0.0041043170 0.006861583 -0.59815890 0.5497339
gene7        0.0056392671 0.006715175  0.83977965 0.4010319
gene8        0.0075179967 0.006381779  1.17804083 0.2387803
gene9       -0.0075242799 0.006423353 -1.17139451 0.2414402
gene10      -0.0028628311 0.006855948 -0.41756897 0.6762623
gene11      -0.0084873719 0.006652622 -1.27579351 0.2020285
gene12       0.0078015271 0.006364037  1.22587699 0.2202450
gene13       0.0042737155 0.006563394  0.65114415 0.5149534
gene14       0.0003330827 0.006167466  0.05400641 0.9569301
gene15      -0.0008043884 0.006419857 -0.12529695 0.9002885
gene16       0.0107164052 0.006445916  1.66251080 0.0964104
gene17       0.0060968062 0.005821314  1.04732472 0.2949498
gene18       0.0098128040 0.005878808  1.66918261 0.0950812
gene19      -0.0064937415 0.006758917 -0.96076652 0.3366696
gene20      -0.0047942275 0.006712493 -0.71422452 0.4750884

• 'Estimate' is the beta coefficient ('beta') - the exponent of this is the odds ratio
• 'Std. Error' is the standard error of the estimate - you will want this to be as low as possible (see https://stats.stackexchange.com/questions/289351/what-does-standard-error-and-prz-mean-in-glm)
• 'z value' is Z value
• 'Pr(>|z|)' is the p-value

[mmh7272]

thank u for ur kindness!