2012-04-08 1 views
0

私は比例データに膨らんだベータ回帰モデルを適合させようとしています。私はパッケージgamlssを使用していて、ファミリーBEINFを指定しています。私は$mu.coefficientsのp値をどのように抽出できるのだろうかと思います。私がfit.3$mu.coefficientsというコマンドを入力したとき(自分のコードの最後に表示されているように)、Mu係数の見積もりしかできませんでした。以下は私のデータの例です。gamlssのmuパラメータのp値

mydata = data.frame(y = c(0.014931087, 0.003880983, 0.006048048, 0.014931087, 
+   0.016469269, 0.013111447, 0.012715517, 0.007981377), index = c(1,1,2,2,3,3,4,4)) 

mydata 
     y  index 
1 0.004517611  1 
2 0.004351405  1 
3 0.007952064  2 
4 0.004517611  2 
5 0.003434018  3 
6 0.003602046  4 
7 0.002370690  4 
8 0.002993016  4 

> library(gamlss) 
> fit.3 = gamlss(y ~ factor(index), family = BEINF, data = mydata) 
> summary(fit.3) 

******************************************************************* 
Family: c("BEINF", "Beta Inflated") 

Call: 
gamlss(formula = y ~ factor(index), family = BEINF, data = mydata) 


Fitting method: RS() 

------------------------------------------------------------------- 
Mu link function: logit 
Mu Coefficients: 

       Estimate Std. Error t value Pr(>|t|) 
(Intercept)  -5.3994  0.1204 -44.858 1.477e-06 
factor(index)2 0.2995  0.1591 1.883 1.329e-01 
factor(index)3 -0.2288  0.1805 -1.267 2.739e-01 
factor(index)4 -0.5017  0.1952 -2.570 6.197e-02 

------------------------------------------------------------------- 
Sigma link function: logit 
Sigma Coefficients: 
      Estimate Std. Error t value Pr(>|t|) 
(Intercept) -4.456  0.2514 -17.72 4.492e-07 
------------------------------------------------------------------- 
Nu link function: log 
Nu Coefficients: 
      Estimate Std. Error t value Pr(>|t|) 
(Intercept) -21.54  10194 -0.002113 0.9984 

------------------------------------------------------------------- 
Tau link function: log 
Tau Coefficients: 
      Estimate Std. Error t value Pr(>|t|) 
(Intercept) -21.63  10666 -0.002028 0.9984 

------------------------------------------------------------------- 
No. of observations in the fit: 8 
Degrees of Freedom for the fit: 7 
     Residual Deg. of Freedom: 1 
         at cycle: 12 

Global Deviance:  -93.08548 
      AIC:  -79.08548 
      SBC:  -78.52938 
******************************************************************* 

fit.3$mu.coefficients 
    (Intercept) factor(index)2 factor(index)3 factor(index)4 
    -5.3994238  0.2994738  -0.2287571  -0.5016511 

本当にありがとうございます。

+2

上記のモデルのために、このようなsummary.gamlssで保存するオプションを使用しますが、それらは、コードと混合されています結果を表示するには: 最も簡単なのはおそらくその関数 を変更して、それが計算するものも返すようにすることです。 –

答えて

1

計算がgamlss ::: summary.gamlss`、 `である

fit.3 = gamlss(y ~ factor(index), family = BEINF, data = mydata)  
sfit.3<-summary(fit.3, save=TRUE) 
sfit.3$mu.coef.table 
sfit.3$sigma.coef.table 
#to get a list of all the slots in the object 
str(sfit.3) 
0
fit.3 = gamlss(y ~ factor(index), family = BEINF, data = mydata.ex)  
sfit.3<-summary(fit.3, save=TRUE) 
fit.3$mu.coefficients 
sfit.3$coef.table # Here use Brackets [] 
estimate.pval<-data.frame(Intercept=sfit.3$coef.table[1,1],pvalue=sfit.3$coef.table[1,4], 
          "factor(index)^2"=sfit.3$coef.table[2,1] ,pvalue=sfit.3$coef.table[2,4], 
          "factor(index)^3"=sfit.3$coef.table[3,1] ,pvalue=sfit.3$coef.table[3,4], 
         "factor(index)^4"=sfit.3$coef.table[4,1] ,pvalue=sfit.3$coef.table[4,4]) 
関連する問題