[R] for help about logistic regression model
Aimin Yan
aiminy at iastate.edu
Tue Nov 21 19:08:57 CET 2006
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.
More information about the R-help
mailing list