作者ching0629 (Syameroke)
看板Statistics
标题Re: [程式] 在R里面用多组资料作线性回归
时间Wed Jan 9 19:53:23 2013
※ 引述《swwf (十字路口)》之铭言：
: ------------------------------------------------------------------------
: [软体程式类别]:R
: [程式问题]: 资料处理、回归
: [软体熟悉度]:
: 低(1~3个月)
: [问题叙述]:
:
: 是这样，我有多组数据需要跑线性回归，y~ ax+b
: 如果只是单纯跑一组的话我还会弄，但如果是很多分组分层的话就不知道怎麽写回圈
: 原始数据是X(时间)，和Y(存活比率)作简单直线回归，
: 现在如果我数据有多项不同的处理，整理结果如下样式
: 处理A 处理B 处理C 重复 时间 存活比率
: 7个层级 2个层级 6个 4个 X1~Xn Y1~ Yn
: 要注意的是，每个组合的时间n是不一样的观察数量，有的组合观察数少(早早全死光)，
: 有的高(等很久才死完)
: 现在想要跑这7*2*6*4个直线回归，把作回归得来的值给输出，
: 并把这些直线回归作图，请问我的回圈程式要怎麽写?
: 搞不清楚要用apply还是tapply...
: 请详尽叙述遭遇到的问题，可能的话，分点叙述你要处理的流程
: [程式范例]:
: 虽然张贴程式很可怕，但基本上有些程式还是要张贴才能解决
: for (i in treat1)
: {for (j in treat2)
: {for (k in treat3)
: {for (l in rep)
: {lm.reg=lm(x~y)
: }
: }
: }
: }
: summary[i,j,k,l](lm.reg)
: plot[i,j,k,l](x,y)
: -----------------------------------------------------------------------------
: 不好意思，我是R和写程式的新手，真的不太知道要怎麽样跑这种回圈，另外想要一页
: 多图把结果的图都画一下，想请教各位大大指导我该如何写这回圈，谢谢。
: 打扰各位不好意思 >_<
问题实在问的很不清楚,我就按照我假设你在想的做吧
你想要做7*2*4*6个回归,而输出的项目是每个回归的"斜率"及"标准误"
如果你还想输出更多东西就自己写在回圈内吧
Output=matrix(NA,nrow=7*2*4*6,ncol=6)
data$treat1=factor(data$treat1)
data$treat2=factor(data$treat2)
data$treat3=factor(data$treat3)
data$rep=factor(data$rep)
for (i in 1:length(levels(treat1))) {
for (j in 1:length(levels(treat2))) {
for (k in 1:length(levels(treat3))) {
for (l in 1:length(levels(rep))) {
dataset=data[data$treat1==levels(treat1)[i]&data$treat2==levels(treat2)[j]&data
$treat3==levels(treat3)[k]&data$rep==levels(rep)[l],]
n=2*4*6*(i-1)+4*6*(j-1)+6*(k-1)+l
Output[n,1]=levels(treat1)[i]
Output[n,2]=levels(treat2)[j]
Output[n,3]=levels(treat3)[k]
Output[n,4]=levels(rep)[l]
if (sum(dataset$x)==0) {
Output[n,5]=NA
Output[n,6]=NA
}
else {
Output[n,5]=lm(y~x,data=dataset)$coefficients[2]
Output[n,6]=summary(lm(y~x,data=dataset))$coefficients[2,2]
}
}
}
}
}
colnames(Output)=c("treat1","treat2","treat3","rep","coefficients","se")
最後Output这项物件就有你要的结果
共有7*2*4*6列
第一行是treat1
第二行是treat2
第三行是treat3
第四行是rep
第五行是符合前4项条件的资料的斜率
第六行是符合前4项条件的资料的斜率的标准误
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 61.64.175.188
※ 编辑: ching0629 来自: 61.64.175.188 (01/09 19:54)
1F：推 swwf:谢谢!!! 01/10 00:14
2F：→ swwf:图的话就无所谓一定要一页全摆完，但想要批次画系列图不会 01/10 00:15
※ 编辑: ching0629 来自: 27.105.2.161 (01/10 00:19)
3F：推 swwf:发现一个问题是，有的数据是存活率直接0，结果就没X值，因为 01/10 00:39
4F：→ swwf:只有一个点，我就没把它放进资料中，因为无存活时间@@ 这样造 01/10 00:40
5F：→ swwf:成各个数据不一样长@@，不知道要怎麽设定跳过那些n/a的 @@ 01/10 00:40
6F：→ ching0629:不太知道你说的是什麽,假设有些列X是0吧,帮你跳过了 01/10 07:59
7F：→ ching0629:如果是NA什麽的,先用x[is.na(x)==1]<-0,把NA设定为0 01/10 08:01
※ 编辑: ching0629 来自: 210.60.122.8 (01/10 08:02)
8F：推 swwf:结果跟我说 'closure' 类型的物件无法具有子集合<--请教何解? 01/10 13:19
9F：→ gsuper:as.factor(X) 01/11 20:14