[R] for help about logistic regression model

[Douglas Bates]

On 11/20/06, Aimin Yan <aiminy at iastate.edu> wrote:
> 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

but apparently that isn't the entire data set.  Can you tell us how
many observations there are in total and how many levels of p and aa
%in% p?  Better yet, could you make the data available, perhaps at a
web site, so we can try fitting the models and see what happens?

I would recommend fitting generalized linear mixed models using the
Laplace approximation to the likelihood rather than PQL.  The Laplace
approximation is now the default method for the lmer function in
package lme4.  The corresponding call would be

library(lme4)
mp5.null <- lmer(Y ~ 1 + (1|p/aa), p5, binomial)

or, to see the progress of the iterations

mp5.null <- lmer(Y ~ 1 + (1|p/aa), p5, binomial, control = list(usePQL
= FALSE, msV = 1))

> 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]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>