Logit in R

Question

I read in the dataset designated Brass and I need to find logit log log (p / 1-p) for 3 countries for each age and speak against the 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 logit in R as

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

but when I try to execute values ​​for Sweden, I get an error. Secondly, how can I build a logistic chart for countries to compare them.

Answer 1

Reading 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 ")) 

Receive packages:

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

Percentage of scaling to proportions:

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

Change the data:

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

Schedule all variables (including Brass_Standard ) by age, with the y axis being converted to a logit scale with linear regression readings:

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

Change the data a little differently:

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

Graph data instead of Brass_Standard , not vs. Age : convert x and y to logit scales and add linear regression 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") 

If you need to get the coefficients of these settings, 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 this might be useful too.

(I hope I will not do my homework for you ...)