I have a dataset like this:

  p  aa 
index        x       y         z      sdx     sdy      sdz   delta     as 
  ms        cur       sc
1 821p MET     1 -5.09688 32.8830 -5.857620 1.478200 1.73998 0.825778 
13.7883 126.91 92.37 -0.1320180 111.0990
2 821p THR     2 -4.07357 28.6881 -4.838430 0.597674 1.37860 1.165780 
13.7207  64.09 50.72 -0.0977129  98.5319
3 821p GLU     3 -5.86733 30.4759 -0.687444 1.313290 1.99912 0.895972 
22.7940  82.73 75.30 -0.0182231  65.0617
4 821p TYR     4 -1.35333 26.9128 -0.491750 1.038350 1.37449 2.285180 
44.7348   7.60 17.17  0.2367850  75.3691
5 821p LYS     5 -4.27189 27.7594  6.272780 1.205650 1.20123 1.633780 
53.3304  41.98 57.68  0.1305950  91.1431
6 821p LEU     6  0.05675 27.5178  6.309750 1.370120 0.64664 1.656920 
27.4681   0.00  0.00  0.0000000  94.0851

here p is random effect, and aa is nested in p
I do like this:
p5 <- read.csv("p_5_angle.csv", header=T, sep=",")
Y<-p5$sc>=90 # probability of pointing inward
library(MASS)

mp5.null <- glmmPQL(Y~1,data=p5,random=~1|p/aa,family=binomial(logit))
summary(mp5.null)

mp5.full<-glmmPQL(Y~as*ms*cur,data=p5,random=~1|p/aa,family=binomial(logit))
summary(mp5.full)

But output give me

Linear mixed-effects model fit by maximum likelihood
  Data: p5
   AIC BIC logLik
    NA  NA     NA

Random effects:
  Formula: ~1 | p
         (Intercept)
StdDev:   0.1165222

  Formula: ~1 | aa %in% p
         (Intercept)  Residual
StdDev:   0.6551498 0.9735705

Variance function:
  Structure: fixed weights
  Formula: ~invwt
Fixed effects: Y ~ 1
                  Value Std.Error  DF   t-value p-value
(Intercept) -0.1256839 0.1117682 938 -1.124505  0.2611

Standardized Within-Group Residuals:
        Min         Q1        Med         Q3        Max
-1.4693363 -0.8816572 -0.5038361  0.9541089  2.0939872

Number of Observations: 1030
Number of Groups:
         p aa %in% p
         5        92

AIC,BIC,LogLik is NA.

Does anybody know why?

thanks,

Aimin Yan




	[[alternative HTML version deleted]]

