4
$\begingroup$

I have 4 groups of animals that did a task on four consecutive days. The variable of interest was latency, which was censored at 60 seconds. I fit a mixed effects linear censored regression model as follows:

library(censReg)
m <- censReg(Latency ~ Group * Day,
  data = data, method = "BHHH", right = 60)
summary(m)

Observations:
         Total  Left-censored     Uncensored Right-censored 
           608              0            482            126 

Coefficients:
             Estimate Std. error t value  Pr(> t)    
(Intercept)  59.64948    4.80290  12.419  < 2e-16 ***
Group2        1.82896    6.22570   0.294  0.76893    
Group3      -26.53572    6.32002  -4.199 2.68e-05 ***
Group4       -4.24361    6.61237  -0.642  0.52102    
Day         -11.98992    1.39992  -8.565  < 2e-16 ***
Group2:Day   -0.63445    1.66085  -0.382  0.70246    
Group3:Day    6.56444    2.05358   3.197  0.00139 ** 
Group4:Day    3.52366    2.42819   1.451  0.14674    
logSigmaMu    2.13480    0.19363  11.025  < 2e-16 ***
logSigmaNu    3.04728    0.03928  77.577  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

BHHH maximisation, 53 iterations
Return code 8: successive function values within relative tolerance limit (reltol)
Log-likelihood: -2302.083 on 10 Df



library(parameters)
model_parameters(m)
         Parameter Coefficient         SE   CI     CI_low    CI_high          t df_error             p
    1  (Intercept)  59.6494822 4.80289885 0.95  50.216882  69.082082 12.4194750      598  1.140841e-31
    2       Group2   1.8289648 6.22570236 0.95 -10.397934  14.055864  0.2937765      598  7.690307e-01
    3       Group3 -26.5357235 6.32002420 0.95 -38.947865 -14.123582 -4.1986743      598  3.093431e-05
    4       Group4  -4.2436141 6.61237005 0.95 -17.229905   8.742677 -0.6417690      598  5.212691e-01
    5          Day -11.9899227 1.39992101 0.95 -14.739282  -9.240563 -8.5647137      598  9.194472e-17
    6   Group2:Day  -0.6344538 1.66084945 0.95  -3.896261   2.627353 -0.3820056      598  7.025929e-01
    7   Group3:Day   6.5644385 2.05358314 0.95   2.531327  10.597550  3.1965779      598  1.464037e-03
    8   Group4:Day   3.5236628 2.42818993 0.95  -1.245154   8.292479  1.4511479      598  1.472630e-01
    9   logSigmaMu   2.1347957 0.19362696 0.95   1.754524   2.515067 11.0253019      598  7.364901e-26
    10  logSigmaNu   3.0472752 0.03928047 0.95   2.970131   3.124420 77.5773541      598 2.539111e-314

I don't know how to do a post-hoc analysis to compare all groups with one another on each day. I always use the emmeans package to get contrasts, but it doesn't work for a censReg object. Do you have any suggestions to get around that? Easy solutions are most welcome :-) You can get the data in a csv file here.


As suggested, I tried with the qdrg() funciton from the emmeans package.

rg <- qdrg(Latency ~ Group * Day, data=data, coef=coef(m),vcov=vcov(m), df=m$df.residual)

'emmGrid' object with variables:
        Group = 1, 2, 3, 4
        Day = 2.5

Although the qdrg() funciton seemed to work, the rg grid did not work with the emmenas function.

emmeans(rg, pairwise ~ Group|Day,  at = list(Day = c(1,2,3, 4)), adjust="none")

Error in X[ii, ii, drop = FALSE] %*% y[ii] : non-conformable arguments

I don't know how to fix it.

$\endgroup$
4
  • 1
    $\begingroup$ Did you try emmeans::qdrg()? $\endgroup$
    – Russ Lenth
    Jul 13, 2022 at 19:45
  • $\begingroup$ @RussLenth Yes, but I couldn't make it work. I was just hitting a wall. I guess there is a way to push it through qdrg(), but I won't be able to do it on my own. $\endgroup$
    – Uki Buki
    Jul 14, 2022 at 7:07
  • 1
    $\begingroup$ I think you need to subset the coefficients and vcov, to exclude the last two elements which are extra parameters. That is, use coef(m)[1:8] and vcov(m)[1:8, 1:8]. Also, you need to put the at spec in theqdrg call, because that is part of the construction of the reference grid. It will be ignored by emmeans. (I thought I had already commented this, not sure what happened). $\endgroup$
    – Russ Lenth
    Jul 14, 2022 at 18:06
  • 1
    $\begingroup$ @Russ Lenth. Indeed. Subsetting coefficients did the trick. Thanks. $\endgroup$
    – Uki Buki
    Jul 17, 2022 at 13:51

1 Answer 1

0
$\begingroup$

To make qdrg() work, I needed to subset the coefficients.

coef <- coef(m)[1:8]
vcov <- vcov(m)[1:8, 1:8]
rg <- qdrg(Latency ~ Group * Day, data=data, coef=coef, vcov=vcov, df=m$df.residual, at = list(Day = c(1,2,3,4)))
emmeans(rg, pairwise ~ Group|Day, adjust="none")
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.