4

我阅读了标记为黄铜的数据集，我需要找到 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 曲线以进行比较。

4

1 回答 1

4

读取数据：

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 轴转换为 logit 标度，显示线性回归拟合：

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

绘制数据 vs.Brass_Standard而不是 vs. Age：将 x 和 y 都转换为 logit 尺度，并再次添加线性回归拟合。

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

（我希望我没有为你做作业……）

于 2012-12-08T01:22:07.010 回答