[R-sig-ME] glmmPQL inquiry

[Sundar Dorai-Raj]

kennedy otwombe said the following on 12/6/2007 12:33 PM:
> Dear R users,
> 
> I have longitudinal data that is all binary and i have run it in SAS using the following code without a problem:
> 
> proc glimmix data=navs;
> class id;
> model y(event='1')=x1 x2 x3/solution distribution=binary;
> random intercept/subject=id type=cs;
> run;
> 
> I have also written a code in R for the same analysis but i am getting the following error message (iteration 1
> Error in switch(mode(x), "NULL" = structure(NULL, class = "formula"),  : 
>   invalid formula). 
> 
> My code in R reads as follows:
> 
>> fit<-glmmPQL(y~x1+x2+x3, random=~1|id, family=binomial, data=navs)
>> summary(fit)
> 
> I am not sure where the problem lies but i realise that glmmPQL does not seem to cater for binary distributions and i aint sure how to model the probability of the event Y=1 which is my interest in the data i am assuming. My data looks as follows:
> 
>   t id y x0 x1 x2 x3
>   1  1 0  0  0  1  1
>   2  1 0  0  1  0  1
>   3  1 0  0  1  0  1
>   1  2 0  0  0  1  0
>   2  2 0  1  1  1  0
>   3  2 0  1  1  1  0
> 
> I will appreciate any ideas from this network.
> 
>  
> Kennedy.N.Otwombe
> School of Statistics & Actuarial Science 
> University of the Witwatersrand 
> Private Bag 3 
> Wits 2050 
> Johannesburg 
> South Africa
> 

Hi, Kennedy,

You have not given us enough information:

1. Result from sessionInfo (includes MASS version assuming glmmPQL is 
from MASS, nlme version, and R version)
2. Realistic data - your example data has y == 0 always. You will never 
produce a realistic model for these observations. And if I fake it and 
assign some ones to "y" the code you provide works.

library(MASS)
navs <- read.table(con <- textConnection("t id y x0 x1 x2 x3
   1  1 0  0  0  1  1
   2  1 0  0  1  0  1
   3  1 0  0  1  0  1
   1  2 0  0  0  1  0
   2  2 0  1  1  1  0
   3  2 1  1  1  1  0"), header = TRUE) ## last row has a 1
close(con)
fit <- glmmPQL(y~x1+x2+x3, random=~1|id, family=binomial, data=navs)

 > summary(fit)
Linear mixed-effects model fit by maximum likelihood
  Data: navs
   AIC BIC logLik
    NA  NA     NA

Random effects:
  Formula: ~1 | id
          (Intercept)  Residual
StdDev: 8.009602e-05 0.5773503

Variance function:
  Structure: fixed weights
  Formula: ~invwt
Fixed effects: y ~ x1 + x2 + x3
                 Value Std.Error DF       t-value p-value
(Intercept) -61.13214   6171227  2 -9.905993e-06       1
x1           30.56607   2631420  2  1.161581e-05       1
x2           30.56607   4160641  2  7.346481e-06       1
x3            0.00000   3721390  0  0.000000e+00     NaN
  Correlation:
    (Intr) x1     x2
x1 -0.853
x2 -0.944  0.632
x3 -0.905  0.707  0.894

Standardized Within-Group Residuals:
           Min            Q1           Med            Q3           Max
-1.732051e+00 -4.715395e-16 -2.827740e-16  4.694641e-17  1.732051e+00

Number of Observations: 6
Number of Groups: 2
Warning message:
In pt(q, df, lower.tail, log.p) : NaNs produced
 > sessionInfo()
R version 2.6.1 (2007-11-26)
i386-pc-mingw32

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United 
States.1252;LC_MONETARY=English_United 
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] nlme_3.1-86 MASS_7.2-38

loaded via a namespace (and not attached):
[1] grid_2.6.1     lattice_0.17-2
 >