glmmPQL inquiry

Hi Sundar/Ken,
I would like to clarify my inquiry. The data structure i gave was just a sample of my data. I actually have 50 subjects each with three binary observations (hence 150 observations). So the error i initially gave emanates from the analysis involving all the 50 subjects. It just turns out that for the sample i gave, the first two subjects have only 0 as their entries but this varies as you look through the data.
I am using the latest version of R i.e 2.6.1 and i downloaded the nlme version currently available on the CRAN website. I am aslo assuming my random part is the intercept.
Hope this clarifies my inquiry.

Kennedy.N.Otwombe
School of Statistics & Actuarial Science 
University of the Witwatersrand 
Private Bag 3 
Wits 2050 
Johannesburg 
South Africa



----- Original Message ----
From: Sundar Dorai-Raj <sundar.dorai-raj at pdf.com>
To: kennedy otwombe <notwombe at yahoo.com>
Cc: R-sig-mixed-models at r-project.org
Sent: Thursday, December 6, 2007 10:57:55 PM
Subject: Re: [R-sig-ME] glmmPQL inquiry



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

____________________________________________________________________________________
Looking for last minute shopping deals?