4

真ちゅうというラベルの付いたデータセットを読み取り、各年齢の3か国のロジット関数ログ（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<-function(x) log(x/(1-x))

しかし、たとえばスウェーデンの値で実行しようとすると、エラーが発生します。次に、各国のロジット曲線をプロットして比較するにはどうすればよいですか。

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 軸をロジット スケールに変換し、線形回帰適合を示します。

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_Standardvs.ではなく vs.データをプロット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も役立つようです。

（私があなたのために宿題をしていないことを願っています...）

于 2012-12-08T01:22:07.010 に答える