[R-sig-ME] glmmPQL inquiry
Sundar Dorai-Raj
sundar.dorai-raj at pdf.com
Thu Dec 6 21:57:55 CET 2007
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
>
More information about the R-sig-mixed-models
mailing list