# [R] random effects model

arun smartpink111 at yahoo.com
Mon Jan 7 20:14:44 CET 2013

```HI,

Regarding question:2) Have you checked summary(m1)?
data(seizure)
## Diggle, Liang, and Zeger (1994) pp166-168, compare Table 8.10
seiz.l <- reshape(seizure,
varying=list(c("base","y1", "y2", "y3", "y4")),
v.names="y", times=0:4, direction="long")
seiz.l <- seiz.l[order(seiz.l\$id, seiz.l\$time),]
seiz.l\$t <- ifelse(seiz.l\$time == 0, 8, 2)
seiz.l\$x <- ifelse(seiz.l\$time == 0, 0, 1)
m1 <- geese(y ~ offset(log(t)) + x + trt + x:trt, id = id,
data=seiz.l, corstr="exch", family=poisson)
summary(m1)
#Call:
#geese(formula = y ~ offset(log(t)) + x + trt + x:trt, id = id,
#   data = seiz.l, family = poisson, corstr = "exch")
#
#Mean Model:
# Variance to Mean Relation: poisson

#Coefficients:
#             estimate    san.se       wald         p
#(Intercept)  1.34760922 0.1573571 73.3423807 0.0000000
#x            0.11183602 0.1159304  0.9306116 0.3347040
#trt          0.02753449 0.2217878  0.0154127 0.9011982
#x:trt       -0.10472579 0.2134448  0.2407334 0.6236769
--------------------------------------------------------------------
#p is the colname with the p values for the Coefficients
summary(m1)\$mean["p"]
#                    p
#(Intercept) 0.0000000
#x           0.3347040
#trt         0.9011982
#x:trt       0.6236769

You didn't say specifically whether you got NA's in example data or your actual data.  I am getting the p-values in R 2.15.

1) I think it should work with binary response variable.
Another example in the same documentation:
data(ohio)
str(ohio)
#'data.frame':    2148 obs. of  4 variables:
# \$ resp : int  0 0 0 0 0 0 0 0 0 0 ...
# \$ id   : int  0 0 0 0 1 1 1 1 2 2 ...
# \$ age  : int  -2 -1 0 1 -2 -1 0 1 -2 -1 ...
# \$ smoke: int  0 0 0 0 0 0 0 0 0 0 ...
It is not even factors

fit <- geese(resp ~ age + smoke + age:smoke, id=id, data=ohio,
family=binomial, corstr="exch", scale.fix=TRUE)
summary(fit)\$mean["p"]
#                    p
#(Intercept) 0.00000000
#age         0.01523698
#smoke       0.09478252
#age:smoke   0.42234200

# also tested after converting to factors

ohio1<-ohio
ohio1\$smoke<-factor(ohio1\$smoke)
ohio1\$age<-factor(ohio1\$age)

str(ohio1)
#'data.frame':    2148 obs. of  4 variables:
# \$ resp : int  0 0 0 0 0 0 0 0 0 0 ...
# \$ id   : int  0 0 0 0 1 1 1 1 2 2 ...
# \$ age  : Factor w/ 4 levels "-2","-1","0",..: 1 2 3 4 1 2 3 4 1 2 ...
# \$ smoke: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
fit1<-geese(resp~age+smoke+age:smoke,id=id,data=ohio1,family=binomial,corstr="exch",scale.fix=TRUE)
summary(fit1)\$mean["p"]
#                      p
#(Intercept)  0.00000000
#age-1        0.60555454
#age0         0.45322698
#age1         0.01187725
#smoke1       0.86262269
#age-1:smoke1 0.17239050
#age0:smoke1  0.32223942
#age1:smoke1  0.36686706

A.K.

----- Original Message -----
From: rex2013 <usha.nathan at gmail.com>
To: r-help at r-project.org
Cc:
Sent: Monday, January 7, 2013 6:15 AM
Subject: Re: [R] random effects model

Hi A.K

Below is the comment I get, not sure why.

BP.sub3 is the stacked data without the missing values.

BP.geese3 <- geese(HiBP~time*MaternalAge,data=BP.sub3,id=CODEA,
family=binomial, corstr="unstructured", na.action=na.omit)Error in
`contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels

Even though age has 3 levels; time has 14 years & 21 years; HIBP is a
binary response outcome.

2) When you mentioned summary(m1)\$mean["p"] what did the p mean? i
used this in one of the gee command, it produced NA as answer?

Many thanks

On Mon, Jan 7, 2013 at 5:26 AM, arun kirshna [via R] <
ml-node+s789695n4654795h72 at n4.nabble.com> wrote:

> Hi,
>
> I am  not very familiar with the geese/geeglm().  Is it from
> library(geepack)?
> "
> Can you tell me if I can use the geese or geeglm function with this data
> eg: : HIBP~ time* Age
> Here age is a factor with 3 levels, time: 2 levels, HIBP = yes/no.
>
> #  CODEA Sex MaternalAge Education Birthplace AggScore IntScore Obese14
> #1     1  NA           3         4          1       NA       NA      NA
> #2     3   2           3         3          1        0        0       0
>  # Overweight14 Overweight21 Obese21 hibp14 hibp21
> #1           NA           NA      NA     NA     NA
> #2            0            1       0      0      0
>
> If I understand your new classification:
> BP.stacknormal<- subset(BP_2b,Obese14==0 & Overweight14==0 & Obese21==0 &
> Overweight21==0)
> BP.stackObese <- subset(BP_2b,(Obese14==1& Overweight14==0 &
> Obese14==1&Overweight14==1)|(Obese14==1&Overweight14==1 & Obese21==1 &
> Overweight21==0)|(Obese14==1&Overweight14==0 & Obese21==0 &
> Overweight21==0)|(Obese14==0 & Overweight14==0 & Obese21==1 &
> Overweight21==0)|(Obese14==0 & Overweight14==0 & Obese21==1 &
> Overweight21==1)|(Obese14==0 & Overweight14==1 & Obese21==1
> &Overweight21==1)|(Obese14==1& Overweight14==1 & Obese21==1&
> Overweight21==1)) #check whether there are more classification that fits to
> #Obese
>  BP.stackOverweight <- subset(BP_2b,(Obese14==0 & Overweight14==1 &
> Obese21==0 & Overweight21==1)|(Obese14==0 &Overweight14==1 & Obese21==0 &
> Overweight21==0)|(Obese14==0 & Overweight14==0 & Obese21==0 &
> Overweight21==1))
> BP.stacknormal\$Categ<-"Normal"
> BP.stackObese\$Categ<-"Obese"
> BP.stackOverweight\$Categ <- "Overweight"
>  BP.newObeseOverweightNormal<-na.omit(rbind(BP.stacknormal,BP.stackObese,BP.stackOverweight))
>
>  nrow(BP.newObeseOverweightNormal)
> #[1] 1581
> BP.stack3 <-
> reshape(BP.newObeseOverweightNormal,idvar="CODEA",timevar="time",sep="_",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21"),c("hibp14","hibp21")),v.names=c("Obese","Overweight","hibp"),direction="long")
>
> library(car)
> BP.stack3\$time<-recode(BP.stack3\$time,"1=14;2=21")
>   #  CODEA Sex MaternalAge Education Birthplace AggScore IntScore  Categ
> time
> #8.1     8   2           4         4          1        0        0 Normal
> 14
> #9.1     9   1           3         6          2        0        0 Normal
> 14
>   #  Obese Overweight hibp
> #8.1     0          0    0
>
> Now, your formula: (HIBP~time*Age), is it MaternalAge?
> If it is, it has three values
> unique(BP.stack3\$MaternalAge)
> #[1] 4 3 5
> and for time (14,21) # If it says that geese/geeglm, contrasts could be
> applied with factors>=2 levels, what is the problem?
> If you take "Categ" variable, it also has 3 levels (Normal, Obese,
> Overweight).
>
>  BP.stack3\$MaternalAge<-factor(BP.stack3\$MaternalAge)
>  BP.stack3\$time<-factor(BP.stack3\$time)
>
> library(geepack)
> # Using one of the example datasets:
> data(seizure)
>      seiz.l <- reshape(seizure,
>                        varying=list(c("base","y1", "y2", "y3", "y4")),
>                        v.names="y", times=0:4, direction="long")
>      seiz.l <- seiz.l[order(seiz.l\$id, seiz.l\$time),]
>      seiz.l\$t <- ifelse(seiz.l\$time == 0, 8, 2)
>      seiz.l\$x <- ifelse(seiz.l\$time == 0, 0, 1)
>      m1 <- geese(y ~ offset(log(t)) + x + trt + x:trt, id = id,
>                  data=seiz.l, corstr="exch", family=poisson)
>      summary(m1)
>
>  summary(m1)\$mean["p"]
> #                    p
> #(Intercept) 0.0000000
> #x           0.3347040
> #trt         0.9011982
> #x:trt       0.6236769
>
>
> #If you need the p-values of the scale
>    summary(m1)\$scale["p"]
>  #                   p
> #(Intercept) 0.0254634
>
> Hope it helps.
>
> A.K.
>
>
>
>
>
>
> ----- Original Message -----
> From: rex2013 <[hidden email]<http://user/SendEmail.jtp?type=node&node=4654795&i=0>>
>
> To: [hidden email] <http://user/SendEmail.jtp?type=node&node=4654795&i=1>
> Cc:
> Sent: Sunday, January 6, 2013 4:55 AM
> Subject: Re: [R] random effects model
>
> Hi A.K
>
> Regarding my question on comparing normal/ obese/overweight with blood
> pressure change, I did finally as per the first suggestion of stacking the
> data and creating a normal category . This only gives me a obese not obese
> 14, but when I did with the wide format hoping to  get  a
> obese14,normal14,overweight 14 Vs hibp 21, i could not complete any of the
> models.
> This time I classified obese=1 & overweight=1 as obese itself.
>
> Can you tell me if I can use the geese or geeglm function with this data
> eg: : HIBP~ time* Age
> Here age is a factor with 3 levels, time: 2 levels, HIBP = yes/no.
>
> It says geese/geeglm: contrast can be applied only with factor with 2 or
> more levels. What is the way to overcome this. Can I manipulate the data
> to
> make it work.
>
> I need to know if the demogrphic variables affect change in blood pressure
> status over time?
>
> How to get the p values with gee model?
>
> Thanks
> On Thu, Jan 3, 2013 at 5:06 AM, arun kirshna [via R] <
> [hidden email] <http://user/SendEmail.jtp?type=node&node=4654795&i=2>>
> wrote:
>
> > HI Rex,
> > If I take a small subset from your whole dataset, and go through your
> > codes:
> >  BP.sub<-BP_2b[410:418,c(1,8:11,13)] #deleted the columns that are not
> > needed
> >  BP.stacknormal<- subset(BP.subnew,Obese14==0 & Overweight14==0)
> > BP.stackObese <- subset(BP.subnew,Obese14==1)
> >  BP.stackOverweight <- subset(BP.subnew,Overweight14==1)
> > BP.stacknormal\$Categ<-"Normal14"
> > BP.stackObese\$Categ<-"Obese14"
> > BP.stackOverweight\$Categ <- "Overweight14"
> >
> BP.newObeseOverweightNormal<-rbind(BP.stacknormal,BP.stackObese,BP.stackOverweight)
>
> >
> >  BP.newObeseOverweightNormal
> > #    CODEA Obese14 Overweight14 Overweight21 Obese21 hibp21        Categ
> > #411   541       0            0            0       0      0     Normal14
> > #415   545       0            0            1       1      1     Normal14
> > #418   549       0            0            1       0      0     Normal14
> > #413   543       1            0            1       1      0      Obese14
> > #417   548       0            1            1       0      0 Overweight14
> > BP.newObeseOverweightNormal\$Categ<-
> > factor(BP.newObeseOverweightNormal\$Categ)
> > BP.stack3 <-
> >
> reshape(BP.newObeseOverweightNormal,idvar="CODEA",timevar="time",sep="_",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21")),v.names=c("Obese","Overweight"),direction="long")
>
> >
> > library(car)
> > BP.stack3\$time<-recode(BP.stack3\$time,"1=14;2=21")
> > BP.stack3 #Here Normal14 gets repeated even at time==21.  Given that you
> > are using the "Categ" and "time" #columns in the analysis, it will give
> > incorrect results.
> > #      CODEA hibp21        Categ time Obese Overweight
> > #541.1   541      0     Normal14   14     0          0
> > #545.1   545      1     Normal14   14     0          0
> > #549.1   549      0     Normal14   14     0          0
> > #543.1   543      0      Obese14   14     1          0
> > #548.1   548      0 Overweight14   14     0          1
> > #541.2   541      0     Normal14   21     0          0
> > #545.2   545      1     Normal14   21     1          1
> > #549.2   549      0     Normal14   21     0          1
> > #543.2   543      0      Obese14   21     1          1
> > #548.2   548      0 Overweight14   21     0          1
> > #Even if I correct the above codes, this will give incorrect
> > results/(error as you shown) because the response variable (hibp21) gets
> > #repeated when you reshape it from wide to long.
> >
> > The correct classification might be:
> >  BP.sub<-BP_2b[410:418,c(1,8:11,13)]
> >
> BP.subnew<-reshape(BP.sub,idvar="CODEA",timevar="time",sep="",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21")),v.names=c("Obese","Overweight"),direction="long")
>
> >
> > BP.subnew\$time<-recode(BP.subnew\$time,"1=14;2=21")
> >  BP.subnew<-na.omit(BP.subnew)
> >
> > BP.subnew\$Categ[BP.subnew\$Overweight==1 & BP.subnew\$time==14 &
> > BP.subnew\$Obese==0]<-"Overweight14"
> > BP.subnew\$Categ[BP.subnew\$Overweight==1 & BP.subnew\$time==21 &
> > BP.subnew\$Obese==0]<-"Overweight21"
> >  BP.subnew\$Categ[BP.subnew\$Obese==1 & BP.subnew\$time==14 &
> > BP.subnew\$Overweight==0]<-"Obese14"
> >  BP.subnew\$Categ[BP.subnew\$Obese==1 & BP.subnew\$time==21 &
> > BP.subnew\$Overweight==0]<-"Obese21"
> >  BP.subnew\$Categ[BP.subnew\$Overweight==1 & BP.subnew\$time==21&
> > BP.subnew\$Obese==1]<-"ObeseOverweight21"
> >  BP.subnew\$Categ[BP.subnew\$Overweight==1 & BP.subnew\$time==14&
> > BP.subnew\$Obese==1]<-"ObeseOverweight14"
> > BP.subnew\$Categ[BP.subnew\$Overweight==0 & BP.subnew\$Obese==0
> > &BP.subnew\$time==14]<-"Normal14"
> >  BP.subnew\$Categ[BP.subnew\$Overweight==0 & BP.subnew\$Obese==0
> > &BP.subnew\$time==21]<-"Normal21"
> >
> > BP.subnew\$Categ<-factor(BP.subnew\$Categ)
> > BP.subnew\$time<-factor(BP.subnew\$time)
> > BP.subnew
> > #      CODEA hibp21 time Obese Overweight             Categ
> > #541.1   541      0   14     0          0          Normal14
> > #543.1   543      0   14     1          0           Obese14
> > #545.1   545      1   14     0          0          Normal14
> > #548.1   548      0   14     0          1      Overweight14
> > #549.1   549      0   14     0          0          Normal14
> > #541.2   541      0   21     0          0          Normal21
> > #543.2   543      0   21     1          1 ObeseOverweight21
> > #545.2   545      1   21     1          1 ObeseOverweight21
> > #548.2   548      0   21     0          1      Overweight21
> > #549.2   549      0   21     0          1      Overweight21
> >
> > #NOw with the whole dataset:
> > BP.sub<-BP_2b[,c(1,8:11,13)] #change here and paste the above lines:
> >     # CODEA hibp21 time Obese Overweight    Categ
> > #3.1      3      0   14     0          0 Normal14
> > #7.1      7      0   14     0          0 Normal14
> > #8.1      8      0   14     0          0 Normal14
> > #9.1      9      0   14     0          0 Normal14
> > #14.1    14      1   14     0          0 Normal14
> > #21.1    21      0   14     0          0 Normal14
> >
> > tail(BP.subnew)
> >   #     CODEA hibp21 time Obese Overweight             Categ
> > #8485.2  8485      0   21     1          1 ObeseOverweight21
> > #8506.2  8506      0   21     0          1      Overweight21
> > #8520.2  8520      0   21     0          0          Normal21
> > #8529.2  8529      1   21     1          1 ObeseOverweight21
> > #8550.2  8550      0   21     1          1 ObeseOverweight21
> > #8554.2  8554      0   21     0          0          Normal21
> >
> > summary(lme.1 <- lme(hibp21~time+Categ+ time*Categ,
> > data=BP.subnew,random=~1|CODEA, na.action=na.omit))
> > #Error in MEEM(object, conLin, control\$niterEM) :
> >   #Singularity in backsolve at level 0, block 1
> > #May be because of the reasons I mentioned above.
> >
> > #YOu didn't mention the library(gee)
> > BP.gee8 <- gee(hibp21~time+Categ+time*Categ,
> > data=BP.subnew,id=CODEA,family=binomial,
> > corstr="exchangeable",na.action=na.omit)
> > #Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
> > #Error in gee(hibp21 ~ time + Categ + time * Categ, data = BP.subnew, id
> =
> > CODEA,  :
> >   #rank-deficient model matrix
> > With your codes, it might have worked, but the results may be inaccurate
> > # After running your whole codes:
> >  BP.gee8 <- gee(hibp21~time+Categ+time*Categ,
> > data=BP.stack3,id=CODEA,family=binomial,
> > corstr="exchangeable",na.action=na.omit)
> > #Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
> > #running glm to get initial regression estimate
> >    #        (Intercept)                   time           CategObese14
> >      #    -2.456607e+01           9.940875e-15           2.087584e-13
> >     # CategOverweight14      time:CategObese14 time:CategOverweight14
> >       #    2.087584e-13          -9.940875e-15          -9.940875e-15
> > #Error in gee(hibp21 ~ time + Categ + time * Categ, data = BP.stack3, id
> =
> > CODEA,  :
> >  # Cgee: error: logistic model for probability has fitted value very
> close
> > to 1.
> > #estimates diverging; iteration terminated.
> >
> > In short, I think it would be better to go with the suggestion in my
> > ObeseOverweight14, ObeseOverweight21 etc) as I showed here.
> >
> > A.K.
> >
> >
> >
> >
> >
> >
> >
> >
> > ------------------------------
> discussion
> > below:
> >
>
> > .
> > NAML<
>
> >
>
>
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/random-effects-model-tp4654346p4654776.html
> Sent from the R help mailing list archive at Nabble.com.
>     [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] <http://user/SendEmail.jtp?type=node&node=4654795&i=3>mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
> ______________________________________________
> [hidden email] <http://user/SendEmail.jtp?type=node&node=4654795&i=4>mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
>
> ------------------------------
> below:
> http://r.789695.n4.nabble.com/random-effects-model-tp4654346p4654795.html
> .
>

--
View this message in context: http://r.789695.n4.nabble.com/random-effects-model-tp4654346p4654826.html
Sent from the R help mailing list archive at Nabble.com.
[[alternative HTML version deleted]]

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help