[R] random effects model
arun
smartpink111 at yahoo.com
Tue Jan 15 14:22:09 CET 2013
Hi,
Check these links:
http://comments.gmane.org/gmane.comp.lang.r.ggplot2/6527
https://groups.google.com/forum/#!msg/ggplot2/nfVjxL0DXnY/5zf50zCeZuMJ
A.K.
________________________________
From: Usha Gurunathan <usha.nathan at gmail.com>
To: arun <smartpink111 at yahoo.com>
Cc: R help <r-help at r-project.org>
Sent: Tuesday, January 15, 2013 6:31 AM
Subject: Re: [R] random effects model
Hi AK
Got an error message with
library(ggplot2) > ggplot(BP.stack1,aes(x=factor(HiBP),fill=Obese))+geom_bar(position="fill") Error in rename(x, .base_to_ggplot, warn_missing = FALSE) : could not find function "revalue" > ggplot(BP.stack1,aes(x=factor(HiBP),fill=Overweight))+geom_bar(position="fill") Error in rename(x, .base_to_ggplot, warn_missing = FALSE) : could not find function "revalue"
I got the dot plot, thanks for that.
I have attached some plots, not sure how to interpret, they had unusual patterns.Is it because of missing data? I tried removing the missing data too. They still appeared the same. Do I need to transform the data?
Thanks in advance.
On Tue, Jan 15, 2013 at 8:54 AM, arun <smartpink111 at yahoo.com> wrote:
HI,
>
>
>BP_2b<-read.csv("BP_2b.csv",sep="\t")
>BP_2bNM<-na.omit(BP_2b)
>
>BP.stack3 <- reshape(BP_2bNM,idvar="CODEA",timevar="time",sep="",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21"),c("hibp14","hibp21")),v.names=c("Obese","Overweight","HiBP"),times=factor(c(1,2)),direction="long")
>library(car)
>BP.stack3$Obese<- recode(BP.stack3$Obese,"1='Obese';0='Not Obese'")
>BP.stack3$Overweight<- recode(BP.stack3$Overweight,"1='Overweight';0='Not Overweight'")
>
>library(ggplot2)
>ggplot(BP.stack3,aes(x=factor(HiBP),fill=Obese))+geom_bar(position="fill")
>ggplot(BP.stack3,aes(x=factor(HiBP),fill=Overweight))+geom_bar(position="fill")
>
>You could try lmer() from lme4.
>library(lme4)
>fm1<-lmer(HiBP~time+(1|CODEA), family=binomial,data=BP.stack3) #check codes, not sure
>print(dotplot(ranef(fm1,post=TRUE),
> scales = list(x = list(relation = "free")))[[1]])
>qmt1<- qqmath(ranef(fm1, postVar=TRUE))
>print(qmt1[[1]])
>
>
>A.K.
>
>
>
>
>
>________________________________
>From: Usha Gurunathan <usha.nathan at gmail.com>
>To: arun <smartpink111 at yahoo.com>
>Cc: R help <r-help at r-project.org>
>Sent: Monday, January 14, 2013 6:32 AM
>
>Subject: Re: [R] random effects model
>
>
>Hi AK
>
>I have been trying to create some plots. All being categorical variables, I am not getting any luck with plots. The few ones that have worked are below:
>
>barchart(~table(HiBP)|Obese,data=BP.sub3) ## BP.sub3 is the stacked data without missing values
>
>barchart(~table(HiBP)|Overweight,data=BP.sub3)
>
>plot(jitter(hibp14,factor=2)~jitter(Obese14,factor=2),col="gray",cex=0.7, data=Copy.of.BP_2) ## Copy.of.BP_2 is the original wide format
>
>## not producing any good plots with mixed models as well.
>summary(lme.3 <- lme(HiBP~time, data=BP.sub3,random=~1|CODEA, na.action=na.omit))
>anova(lme.3)
>head(ranef(lme.3))
>print(plot(ranef(lme.3))) ##
>
>Thanks for any help.
>
>
>
>
>
>On Mon, Jan 14, 2013 at 4:33 AM, arun <smartpink111 at yahoo.com> wrote:
>
>
>>
>>
>>HI,
>>
>>I think I mentioned to you before that when you reshape the
>>columns excluding the response variable, response variable gets repeated
>>(in this case hibp14 or hibp21) and creates the error"
>>
>>
>>I run your code, there are obvious problems in the code so I didn't reach up to BP.gee
>>
>>
>>BP_2b<-read.csv("BP_2b.csv",sep="\t")
>>BP.stack3 <- reshape(BP_2b,idvar="CODEA",timevar="time",sep="_",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21")),v.names=c("Obese","Overweight"),times=factor(c(1,2)),direction="long")
>>
>>
>>BP.stack3 <- transform(BP.stack3,CODEA=factor(CODEA),Sex=factor(Sex,labels=c("Male","Female")),MaternalAge=factor(MaternalAge,labels=c("39years or less","40-49 years","50 years or older")),Education=factor(Education,labels=c("Primary/special","Started secondary","Completed grade10", "Completed grade12", "College","University")),Birthplace=factor(Birthplace,labels=c("Australia","Other English-speaking","Other")))
>>
>> BP.stack3$Sex <- factor(BP.stack3$Sex,levels=levels(BP.stack3$Sex)[c(2,1)])
>>
>> BP.sub3a <- subset(BP.stack3,subset=!(is.na(Sex)| is.na(Education)|is.na(Birthplace)|is.na(Education)|is.na(hibp14)| is.na(hibp21)))
>> nrow(BP.sub3a)
>>#[1] 3364
>> BP.sub5a <- BP.sub3a[order(BP.sub3a$CODEA),] # your code was BP.sub5a <- BP.sub3a[order(BP.sub5a$CODEA),]
>> ^^^^^ was not defined before
>>#Next line
>>BPsub3$Categ[BPsub6$Overweight==1&BPsub3$time==1&BPsub3$Obese==0]<- "Overweight14" #It should be BP.sub3 and what is BPsub6, it was not defined previously.
>>#Error in BPsub3$Categ[BPsub6$Overweight == 1 & BPsub3$time == 1 & BPsub3$Obese == :
>> #object 'BPsub3' not found
>>
>>
>>
>>
>>
>>
>>
>>A.K.
>>
>>
>>________________________________
>>From: Usha Gurunathan <usha.nathan at gmail.com>
>>To: arun <smartpink111 at yahoo.com>
>>
>>Sent: Sunday, January 13, 2013 1:51 AM
>>
>>Subject: Re: [R] random effects model
>>
>>
>>
>>HI AK
>>
>>Thanks a lot for explaining that.
>>
>>1. With the chi sq. ( in order to find out if the diffce is significant between groups) do I have create a separate excel file and make a dataframe.How do I go about it?
>>
>>I have resent a mail to Jun Yan at a difft email ad( first add.didn't work, mail not delivered).
>>
>>2. With my previous query ( reg. Obese/Overweight/ Normal at age 14 Vs change of blood pressure status at 21), even though I had compromised without the age-specific regression, but I am still keen to explore why the age-specific regression didn't work despite it looking okay. I have given below the syntax. If you get time, could you kindly look at it and see if it could work by any chance? Apologies for persisting with this query.
>>
>>
>>BP.stack3 <-
>>reshape(Copy.of.BP_2,idvar="CODEA",timevar="time",sep="_",varying=list(c("Obese14","Obese21"),c("Overweight14","Overweight21")),v.names=c("Obese","Overweight"),times=factor(c(1,2)),direction="long
>>BP.stack3
>>head(BP.stack3)
>>tail(BP.stack3)
>>
>> names(BP.stack3)[c(2,3,4,5,6,7)] <-
>>c("Sex","MaternalAge","Education","Birthplace","AggScore","IntScore")
>>
>>BP.stack3 <-
>>transform(BP.stack3,CODEA=factor(CODEA),Sex=factor(Sex,labels=c("Male","Female")),MaternalAge=factor(MaternalAge,labels=c("39years
>>or less","40-49 years","50 years or
>>older")),Education=factor(Education,labels=c("Primary/special","Started
>>secondary","Completed grade10", "Completed grade12",
>>"College","University")),Birthplace=factor(Birthplace,labels=c("Australia","Other
>>English-speaking","Other")))
>>
>>table(BP.stack3$Sex)
>>BP.stack3$Sex <-
>>factor(BP.stack3$Sex,levels=levels(BP.stack3$Sex)[c(2,1)])
>>
>>levels(BP.stack3$Sex)
>>BP.sub3a <- subset(BP.stack3,subset=!(is.na(Sex)| is.na(Education)|is.na(Birthplace)|is.na(Education)|is.na(hibp14)| is.na(hibp21)))
>>summary(BP.sub3a)
>>BP.sub5a <- BP.sub3a[order(BP.sub5a$CODEA),]
>> BPsub3$Categ[BPsub6$Overweight==1&BPsub3$time==1&BPsub3$Obese==0]
>><- "Overweight14"
>>BPsub3$Categ[BPsub6$Overweight==1&BPsub3$time==2&BPsub3$Obese==0]
>><- "Overweight21"
>>BPsub3$Categ[BPsub3$Obese==1&BPsub3$time==1&BPsub3$Overweight==0|BPsub3$Obese==1&BPsub3$time==1&BPsub3$Overweight==1
>>] <- "Obese14"
>>BPsub3$Categ[BPsub3$Obese==0&BPsub3$time==1&BPsub3 BPsub3$Categ[BPsub6$Overweight==1&BPsub3$time==1&BPsub3$Obese==0]
>><- "Overweight14"$Overweight==0]
>>
>><- "Normal14"
>>BPsub3$Categ[BPsub3$Obese==0&BPsub3$time==2&BPsub3$Overweight==0]
>><- "Normal21"
>>BPsub3$Categ[BPsub3$Obese==1&BPsub3$time==2&BPsub3$Overweight==0|BPsub3$Obese==1&BPsub3$time==2&BPsub3$Overweight==1]
>><- "Obese21"
>>
>>
>>
>>BPsub3$Categ <- factor(BPsub3$Categ)
>>BPsub3$time <- factor(BPsub3$time)
>>summary(BPsub3$Categ)
>>BPsub7 <- subset(BPsub6,subset=!(is.na(Categ)))
>>BPsub7$time <-
>>recode(BPsub7$time,"1=14;2=21")
>>BPsub7$hibp14 <- factor(BPsub7$hibp14)
>>levels(BPsub7$hibp14)
>>levels(BPsub7$Categ)
>>names(BPsub7)
>>head(BPsub7) ### this was looking quite okay.
>>
>>tail(BPsub7)
>>str(BPsub7)
>>
>>library(gee)
>>
>>BP.gee <- geese(hibp14~ time*Categ,
>>data=BPsub7,id=CODEA,family=binomial,
>>corstr="exchangeable",na.action=na.omit)
>>
>>Thanks again.
>>
>>
>>
>>On Sun, Jan 13, 2013 at 1:22 PM, arun <smartpink111 at yahoo.com> wrote:
>>
>>HI,
>>>
>>>table(BP_2b$Sex) #original dataset
>>># 1 2
>>>#3232 3028
>>> nrow(BP_2b)
>>>#[1] 6898
>>> nrow(BP_2bSexNoMV)
>>>#[1] 6260
>>> 6898-6260
>>>#[1] 638 #these rows were removed from the BP_2b to create BP_2bSexNoMV
>>>BP_2bSexMale<-BP_2bSexNoMV[BP_2bSexNoMV$Sex=="Male",]
>>> nrow(BP_2bSexMale)
>>>#[1] 3232
>>> nrow(BP_2bSexMale[!complete.cases(BP_2bSexMale),]) #Missing rows with Male
>>>#[1] 2407
>>> nrow(BP_2bSexMale[complete.cases(BP_2bSexMale),]) #Non missing rows with Male
>>>#[1] 825
>>>
>>>
>>>You did the chisquare test on the new dataset with 6260 rows, right.
>>>I removed those 638 rows because these doesn't belong to either male or female, but you want the % of missing value per male or female. So, I thought this will bias the results. If you want to include the missing values, you could do it, but I don't know where you would put that missing values as it cannot be classified as belonging specifically to males or females. I hope you understand it.
>>>
>>>Sometimes, the maintainer's respond a bit slow. You have to sent an email reminding him again.
>>>
>>>Regarding the vmv package, you could email Waqas Ahmed Malik (malik at math.uni-augsburg.de) regarding options for changing the title and the the font etc.
>>>You could also use this link (http://www.r-bloggers.com/visualizing-missing-data-2/ ) to plot missing value (?plot.missing()). I never used that package, but you could try. Looks like it gives more information.
>>>
>>>A.K.
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>________________________________
>>>From: Usha Gurunathan <usha.nathan at gmail.com>
>>>To: arun <smartpink111 at yahoo.com>
>>>Sent: Saturday, January 12, 2013 9:05 PM
>>>
>>>Subject: Re: [R] random effects model
>>>
>>>
>>>Hi A.K
>>>
>>>So it is number of females missing/total female participants enrolled: 72.65%
>>>Number of females missing/total (of males+ females) participants enrolled : 35.14%
>>>
>>>The total no. with the master data: Males: 3232, females: 3028 ( I got this before removing any missing values)
>>>
>>>with table(Copy.of.BP_2$ Sex) ## BP
>>>
>>>
>>>If I were to write a table ( and do a chi sq. later),
>>>
>>>as Gender Study Non study/missing Total
>>> Male 825 (25.53%) 2407 (74.47%) 3232 (100%)
>>> Female 828 (27.35%) 2200 (72.65%) 3028 ( 100%)
>>> Total 1653 4607 6260
>>>
>>>
>>>The problem is when I did
>>>>colSums(is.na(Copy.of.BP_2), the sex category showed N=638.
>>>
>>>I cannot understand the discrepancy.Also, when you have mentioned to remove NA, is that not a missing value that needs to be included in the total number missing. I am a bit confused. Can you help?
>>>
>>>## I tried sending email to gee pack maintainer at the ID with R site, mail didn't go through??
>>>
>>>Many thanks
>>>
>>>
>>>
>>>
>>>
>>>
>>>On Sun, Jan 13, 2013 at 9:17 AM, arun <smartpink111 at yahoo.com> wrote:
>>>
>>>Hi,
>>>>Yes, you are right. 72.655222% was those missing among females. 35.14377% of values in females are missing from among the whole dataset (combined total of Males+Females data after removing the NAs from the variable "Sex").
>>>>
>>>>A.K.
>>>>
>>>>
>>>>
>>>>________________________________
>>>>From: Usha Gurunathan <usha.nathan at gmail.com>
>>>>To: arun <smartpink111 at yahoo.com>
>>>>Cc: R help <r-help at r-project.org>
>>>>Sent: Saturday, January 12, 2013 5:59 PM
>>>>
>>>>Subject: Re: [R] random effects model
>>>>
>>>>
>>>>
>>>>Hi AK
>>>>That works. I was trying to get similar results from any other package. Being a beginner, I was not sure how to modify the syntax to get my output.
>>>>
>>>>lapply(split(BP_2bSexNoMV,BP_
>>>>2bSexNoMV$Sex),function(x) (nrow(x[!complete.cases(x[,-2]),])/nrow(x))*100) #gives the percentage of rows of missing #values from the overall rows for Males and Females
>>>>#$Female
>>>>#[1] 72.65522
>>>>#
>>>>#$Male
>>>>#[1] 74.47401
>>>>
>>>>#iF you want the percentage from the total number rows in Males and Females (without NA's in the the Sex column)
>>>> lapply(split(BP_2bSexNoMV,BP_2bSexNoMV$Sex),function(x) (nrow(x[!complete.cases(x[,-2]),])/nrow(BP_2bSexNoMV))*100)
>>>>#$Female
>>>>#[1] 35.14377
>>>>#
>>>>#$Male
>>>>#[1] 38.45048
>>>>
>>>>How do I interpret the above 2 difft results? 72.66% of values were missing among female participants?? Can you pl. clarify.
>>>>
>>>>Many thanks.
>>>>
>>>>
>>>>On Sun, Jan 13, 2013 at 3:28 AM, arun <smartpink111 at yahoo.com> wrote:
>>>>
>>>>lapply(split(BP_2bSexNoMV,BP_2bSexNoMV$Sex),function(x) (nrow(x[!complete.cases(x[,-2]),])/nrow(x))*100) #gives the percentage of rows of missing #values from the overall rows for Males and Females
>>>>>#$Female
>>>>>#[1] 72.65522
>>>>>#
>>>>>#$Male
>>>>>#[1] 74.47401
>>>>>
>>>>>#iF you want the percentage from the total number rows in Males and Females (without NA's in the the Sex column)
>>>>> lapply(split(BP_2bSexNoMV,BP_2bSexNoMV$Sex),function(x) (nrow(x[!complete.cases(x[,-2]),])/nrow(BP_2bSexNoMV))*100)
>>>>>#$Female
>>>>>#[1] 35.14377
>>>>>#
>>>>>#$Male
>>>>>#[1] 38.45048
>>>>
>>>
>>
>
More information about the R-help
mailing list