1

I am fitting a mixture of beta regressions using the betamix package in R. If I use the example given in the package with code:

data("ReadingSkills", package = "betareg")
set.seed(4040)
rs_mix <- betamix(accuracy ~ iq, data = ReadingSkills, k = 3,
              nstart = 10, extra_components = extraComponent(type = "uniform",
                                                             coef = 0.99, delta = 0.01))
summary(rs_mix)

Running the summary of the fitted betamix object gives the results:

> summary(rs_mix)
$Comp.1
$Comp.1$mean
            Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  1.40342    0.26332  5.3296  9.84e-08 ***
iq           0.82502    0.21630  3.8142 0.0001366 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$Comp.1$precision
            Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  2.68509    0.45435  5.9097 3.427e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


$Comp.2
$Comp.2$mean
             Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  0.502523   0.082476  6.0930 1.108e-09 ***
iq          -0.048415   0.112923 -0.4287    0.6681    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$Comp.2$precision
            Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  4.25160    0.74737  5.6888 1.279e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The Component 2 results show the variable iq is not significant. Is there a way to remove this variable from the summary result? I have tried using summary(rs_mix)$Comp.1 and it gives me the error:

Error in summary(rs_mix)$Comp.1 : 
  $ operator not defined for this S4 class
Divi
  • 1,614
  • 13
  • 23
  • you want to remove it from the summary but keep it in the model? – marbel Feb 24 '16 at 23:08
  • Yes, as `iq` is significant for Component 1, it will remain in the model. – Divi Feb 24 '16 at 23:09
  • Would it work to keep the variables from the precision object, and the subset the mean object, for each component? – marbel Feb 24 '16 at 23:14
  • I am not sure I follow exactly, but only analysing the mean object might work for this example as it fits only intercepts to the precision object. – Divi Feb 24 '16 at 23:17

1 Answers1

2

You can access S4 objects using the @ symbol. It works similar to $ for S3 objects.

smr = summary(rs_mix)
comps = smr@components[[1]]
significance_level = 0.05

lapply(comps, function(x) {
  bool = x[['mean']][,4] < significance_level
  x[['mean']][bool, ]
})

$Comp.1
             Estimate Std. Error  z value     Pr(>|z|)
(Intercept) 1.4034184  0.2633230 5.329647 9.840410e-08
iq          0.8250193  0.2163036 3.814172 1.366404e-04

$Comp.2
    Estimate   Std. Error      z value     Pr(>|z|) 
5.025226e-01 8.247566e-02 6.092981e+00 1.108272e-09 
marbel
  • 7,560
  • 6
  • 49
  • 68