CODEWITHSUNDEEP
rmath
mark
mark

Reputation:

logit function in R

I read in the data set labeled Brass and I need to find the logit function log (p/1-p) for the 3 countries for each age and plot against brass standard. 

dat <- structure(list(Age=c(1L,5L,10L,20L,30L),Brass_Standard=c(85,76.9,75,71.3,65.2),Sweden=c(98.7,98.4,98.2,97.9,97.4),Italy=c(84.8,73.9,72.1,69.9,64.1),Japan=c(96.4,95.2,94.7,93.8,91.7)),.Names=c("Age","Brass_Standard","Sweden","Italy","Japan"),class="data.frame",row.names=c("1","2","3","4","5"))

     Age   Brass_Standard  Sweden  Italy  Japan
1      1             85.0    98.7   84.8   96.4
2      5             76.9    98.4   73.9   95.2
3     10             75.0    98.2   72.1   94.7
4     20             71.3    97.9   69.9   93.8
5     30             65.2    97.4   64.1   91.7

I defined the logit in R as

logit<-function(x) log(x/(1-x))

but when I try to perform on the values for say Sweden, I get an error. Secondly,how do I plot the logit curve for the countries to compare them.

Upvotes: 4

Views: 10261

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226247

Read data:

dat <- read.table(textConnection(
"Age Brass_Standard    Sweden  Italy   Japan
1   1   85.0   98.7   84.8   96.4
2   5   76.9   98.4   73.9   95.2
3   10  75.0   98.2   72.1   94.7
4   20  71.3   97.9   69.9   93.8
5   30  65.2   97.4   64.1   91.7
"))

Get packages:

library(ggplot2)
library(scales)
library(reshape2)

Rescale percentages to proportions:

dat[,-1] <- dat[,-1]/100

Reshape data:

mdat <- melt(dat,id.var="Age")

Plot all variables (including Brass_Standard) vs age, with y-axis transformed to a logit scale, with linear regression fits shown:

qplot(Age,value,data=mdat,colour=variable)+
    scale_y_continuous(trans=logit_trans())+
    geom_smooth(method="lm")+theme_bw()
ggsave("logitplot1.png")

enter image description here

Reshape data slightly differently:

mdat2 <- melt(dat,id.var=c("Age","Brass_Standard"))

Plot data vs. Brass_Standard rather than vs. Age: transform both x and y to logit scales, and add linear regression fits again.

qplot(Brass_Standard,value,data=mdat2,colour=variable)+
    scale_y_continuous(trans=logit_trans())+
    scale_x_continuous(trans=logit_trans())+
    geom_smooth(method="lm")+
    theme_bw()
ggsave("logitplot2.png")

enter image description here

If you need to get the coefficients of these fits I would suggest something like:

library(nlme)
pdat <- with(mdat2,data.frame(Age,variable,
                              logit_Brass_Standard=plogis(Brass_Standard),
                              logit_value=plogis(value)))

fit1 <- lmList(logit_Brass_Standard~logit_value|variable,data=pdat)
coef(fit1)

http://www.demog.berkeley.edu/~eddieh/toolbox.html#BrassMortality looks like it might be useful too.

(I hope I'm not doing your homework for you ...)

Upvotes: 4

Related Questions