r - R中的logit函数

我读了标有黄铜的数据集，我需要找到每个年龄段这3个国家的logit函数日志（p / 1-p），并根据黄铜标准进行绘制。

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

我将R中的logit定义为

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

但是，当我尝试执行瑞典的值时，出现错误。其次，我该如何绘制Logit曲线以供各国比较。

最佳答案

读取数据：

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

获取软件包：

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

将百分比调整为比例：

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

重塑数据：

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

绘制所有变量（包括Brass_Standard）与年龄的关系图，将y轴转换为对数刻度，并显示线性回归拟合：

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

重塑数据略有不同：

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

绘制数据与Brass_Standard而不是与Age的关系：将x和y都转换为对数刻度，并再次添加线性回归拟合。

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

如果您需要获得这些拟合的系数，我会建议类似以下内容：

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看起来也可能有用。

（希望我不要为你做功课...）