[R] for help about logistic regression model

[Aimin Yan]

thanks, Here is data under this link with file name as p_5_angle.csv
http://www.public.iastate.edu/~aiminy/data/
p is protein names(5 proteins)
aa are nested in p(up to 19 levels for each p, some p doesn't have 19 levels)
index is position of aa.
there are only one observation for each position of each aa within p.

consider p as random effect,
since aa is nested in p, so aa is also random effect.

p and aa are qualitative predictors.
x,y,z,sdx,sdy,sdz,delta,as,ms,cur are quantitative predictors.
sc is binary responsible variable(>=90 and <90)

we want to know the effect of p,aa,x,y,z,sdx,sdy,sdz,delta,as,ms,cur) on 
P(sc>=90).

So I consider to use logistic regression model with p and aa as random effect.

Firstly I try to use p,aa,x,y,z,sdx,sdy,sdz,delta,as,ms,cur as predictors, 
but it seems it has too many predictors.
so I use p,aa,as,ms,cur as predictors, but it still doesn't work.

I hope that I could get your help.

thanks,

At 10:58 AM 11/21/2006, you wrote:
>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.