4

パッケージを使用して、データセットに対して正確なロジスティック回帰を実行しましたelrm

通常のロジスティック回帰と比較しています。

通常のロジスティック回帰でブートストラップを実行することができました。取得した関心のある統計は、推定係数と p 値でした。

ただし、出力から必要な係数を引き出すことができないため、elrm ブートストラップを実行できません。

私のデータを使用すると、要約が印刷されます。

Results:
   estimate p-value p-value_se mc_size
M  0.15116 0.06594    0.00443   49000

95% Confidence Intervals for Parameters

     lower     upper
M 0.00156155 0.3647232

ブートストラップを実行するときにこれらの統計を取得できるように、M 推定値と p 値を抽出したいと考えています。さまざまな組み合わせを試して値を取得しようとしましたが、機能しません。

summary(model)$coefficient 
summary(model)$Results
summary(model)$estimate

これらはすべて、要約を再び吐き出すだけです。

elrm の概要から抽出できるかどうかは誰にもわかりませんか?

どんな助けでも大歓迎です。

4

2 に答える 2

3 
library(elrm)
data(drugDat)
drug.elrm <- elrm(formula=recovered/n~sex+treatment,interest=~sex+treatment,r=4, 
    iter=100000,burnIn=1000,dataset=drugDat)

> summary.elrm(drug.elrm)

Call:
[[1]]
elrm(formula = recovered/n ~ sex + treatment, interest = ~sex + 
    treatment, r = 4, iter = 1e+05, dataset = drugDat, burnIn = 1000)


Results:
          estimate p-value p-value_se mc_size
joint           NA 0.14886    0.00173   99000
sex        0.27092 0.69385    0.01204    2649
treatment  0.76739 0.07226    0.00314   13160

95% Confidence Intervals for Parameters

               lower    upper
sex       -0.6217756 1.212499
treatment -0.1216884 1.852346

# If you look at the summary function, it simply outputs formatted results to
# the screen. So instead, we can just work with the original drug.elrm object

names(drug.elrm)
# shows you everything in this object

# to see the p-values
drug.elrm$p.values.se
      joint         sex   treatment 
0.001734482 0.012039701 0.003143006 

# to get the p-value for joint
drug.elrm$p.values.se[[1]]

# now for the CI
drug.elrm$coeffs.ci
               lower    upper
sex       -0.6217756 1.212499
treatment -0.1216884 1.852346

> drug.elrm$coeffs.ci[[1]]
[1] -0.6217756 -0.1216884
> drug.elrm$coeffs.ci[[1]][1]
[1] -0.6217756
>
于 2012-06-29T19:07:51.267 に答える
2

p 値と信頼区間は、モデル オブジェクト自体の一部として保存されているようです。このsummary.elrm関数は、その情報の一部を抽出してフォーマットするだけです。これを理解するために使用した手順を以下に示しますが、要約バージョンは、object$p.valuesおよびobject$coeffs.ciそれぞれのモデル オブジェクト自体にインデックスを付けることです。

#Example from help page
data(utiDat) 
uti.elrm <- elrm(uti/n~age+current+dia+oc+pastyr+vi+vic+vicl+vis,
                 interest=~dia,r=4,iter=5000,burnIn=1000,dataset=utiDat)
#Look at summary
summary(uti.elrm)
#Let's examine the structure of the summary object to see what's in there, i.e.
#what can we extract?
str(summary(uti.elrm))
#Hmm, doesn't look like anything of interest. Let's look at the source code itself
summary.elrm
#Looks like the p.values are stored in the actual model object iself and the summary function
#just formats them for us. The relevant part of the summary code for the p-value is:
#-----
  #inferences = as.data.frame(cbind(round(as.numeric(object$coeffs),5), 
  #                                 round(as.numeric(object$p.values), 5), 
  #                                 round(as.numeric(object$p.values.se),5),
  #                                 object$mc.size)
  #                           )
  #results = data.frame(row.names = names(object$coeffs), inferences)
  #names(results) = c("estimate", "p-value", "p-value_se", "mc_size")
#So, it looks like we can grab the p.values directly from "object"
> uti.elrm$p.values
dia 
0.02225
#And the confidence intervals are also in the object, located here in the summary code:
#-----
  #cat(object$ci.level, "% Confidence Intervals for Parameters\n", 
  #    sep = "")
  #cat("\n")
  #print(object$coeffs.ci)
#So we can extract them thusly:
> uti.elrm$coeffs.ci
lower upper
dia 0.1256211   Inf
于 2012-06-29T19:14:38.083 に答える