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]]
I have a dataset that has these variables: p,aa,x,y, z,sdx,sdy,sdz,delta,as,ms,cur,sc here p is random effect, and aa is nested in p p and aa are both qualitative predictor,all other are quantitative predictor. sc is binary responsible variable. I want to know the effect of p,aa,x,y,z,sdx,sdy,sdz,delta,as,ms,cur on sc. How can fit these predictors into a logistic regression model in R? thanks, Aimin Yan
help.search("logistic")
gives some 20 hits of which at least 3 or so are immediately relevant such
as glm, lrm and polr
Hth, Ingmar
> From: Aimin Yan <aiminy at iastate.edu>
> Date: Mon, 20 Nov 2006 23:57:44 -0600
> To: <r-help at stat.math.ethz.ch>
> Subject: [R] for help about logistic regression model
>
> I have a dataset that has these variables:
>
> p,aa,x,y, z,sdx,sdy,sdz,delta,as,ms,cur,sc
>
>
> here p is random effect, and aa is nested in p
> p and aa are both qualitative predictor,all other are quantitative
> predictor. sc is binary responsible variable.
> I want to know the effect of p,aa,x,y,z,sdx,sdy,sdz,delta,as,ms,cur on sc.
>
> How can fit these predictors into a logistic regression model in R?
>
> thanks,
>
> Aimin Yan
>
> ______________________________________________
> 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.
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.0851but 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. >
I want to do model selection,here is null model and full model. p5.lgm.null<-lmer(Y~1+(1|p/aa),data=p5,family=binomial,control=list(usePQL=FALSE,msV=1)) summary(p5.lgm.null) p5.lgm.full <- lmer(Y ~as*ms*cur+(1|p/aa),data=p5,family=binomial,control=list(usePQL=FALSE,msV=1)) summary(p5.lgm.full) I do like this. > p5.lgm.bwd <- stepAIC(p5.lgm.full, direction="backward") Error in terms.default(object) : no terms component > summary(p5.lgm.bwd) Error: object "p5.lgm.bwd" not found Error in summary(p5.lgm.bwd) : error in evaluating the argument 'object' in selecting a method for function 'summary' > it seems stepAIC doesn't work for lmer model. Is this true? Aimin Yan > p5.lgm.9 <- lmer(Y ~p*aa*index*x*y*z*sdx*sdy*sdz*delta*as*ms*cur+(1|p/aa),data=p5,family=binomial,control=list(usePQL=FALSE,msV=1)) Error: cannot allocate vector of size 1565600 Kb In addition: Warning messages: 1: Reached total allocation of 494Mb: see help(memory.size) 2: Reached total allocation of 494Mb: see help(memory.size) At 12:45 PM 11/21/2006, you wrote:>On 11/21/06, Aimin Yan <aiminy@iastate.edu> wrote: >>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. > >Here are the fits for two of your models using lmer. > >>library(lme4) >Loading required package: Matrix >Loading required package: lattice >>p5 <- read.csv("http://www.public.iastate.edu/~aiminy/data/p_5_angle.csv") >>p5$Y <- p5$sc >= 90 >>(mp5.NULL <- lmer(Y ~ 1|p/aa, p5, binomial, control = list(usePQL = FALSE))) >Generalized linear mixed model fit using Laplace >Formula: Y ~ 1 | p/aa > Data: p5 >Family: binomial(logit link) > AIC BIC logLik deviance >1390 1405 -692 1384 >Random effects: >Groups Name Variance Std.Dev. >aa:p (Intercept) 0.447654 0.66907 >p (Intercept) 0.015078 0.12279 >number of obs: 1030, groups: aa:p, 92; p, 5 > >Estimated scale (compare to 1 ) 0.9736361 > >Fixed effects: > Estimate Std. Error z value Pr(>|z|) >(Intercept) -0.1325 0.1151 -1.151 0.25 >>(mp5.full <- lmer(Y ~ as*ms*cur + (1|p/aa), p5, binomial, control = >>list(usePQL = FALSE))) >Generalized linear mixed model fit using Laplace >Formula: Y ~ as * ms * cur + (1 | p/aa) > Data: p5 >Family: binomial(logit link) > AIC BIC logLik deviance >1278 1327 -628.8 1258 >Random effects: >Groups Name Variance Std.Dev. >aa:p (Intercept) 0.085104 0.29173 >p (Intercept) 0.026769 0.16361 >number of obs: 1030, groups: aa:p, 92; p, 5 > >Estimated scale (compare to 1 ) 0.9833564 > >Fixed effects: > Estimate Std. Error z value Pr(>|z|) >(Intercept) 5.506e-01 1.714e-01 3.213 0.00131 ** >as -3.964e-02 2.322e-02 -1.707 0.08778 . >ms 1.879e-02 2.149e-02 0.874 0.38206 >cur 3.413e-01 6.706e-01 0.509 0.61078 >as:ms 1.091e-04 7.615e-05 1.432 0.15201 >as:cur 8.315e-02 7.069e-02 1.176 0.23951 >ms:cur -4.880e-02 5.372e-02 -0.908 0.36366 >as:ms:cur -3.998e-04 6.602e-04 -0.606 0.54476 >--- >Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > >Correlation of Fixed Effects: > (Intr) as ms cur as:ms as:cur ms:cur >as -0.028 >ms -0.144 -0.960 >cur -0.391 0.070 0.019 >as:ms 0.290 -0.655 0.443 -0.155 >as:cur -0.094 0.401 -0.288 0.211 -0.530 >ms:cur 0.118 0.601 -0.672 -0.473 -0.185 -0.263 >as:ms:cur 0.110 -0.354 0.232 0.001 0.614 -0.874 0.036[[alternative HTML version deleted]]