一位读者问到怎么计算几何均值的置信区间。
我第一想到的就是proc means,但是这只能计算均值的置信区间。所以应该需要进行一些变换。
一:基本定义
首先几何均数的计算方法有以下两种(观测没有重复值,有重复值需要用到加权法)
注意第二个计算公式,因为我们在SAS一般是用log计算,但是SAS中直接用log是以e底,然后后面我们需要用到exp函数,这个也是以e为底;但是上面的公式是lg,这个是以10为底。
所以一开始这让我有些困惑,log10和loge算出来的结果肯定是不一样的。通过搜索相关资料,如下面的图片说到,运用下面的公式,只要保证对数和反对数用的底保持一致就行了(注意下面公式左边还是log的形式,所以看起来和上面的公式不太一样)
困惑解开了,接下来直接根据公式计算就好了。
二:SAS实现
以SASHELP里的class数据集为例,假设我们要对height求几何均值及其对应的置信区间。
先对height求log,然后求出均值,最后对平均的对数值取指数,就得到几何均值,然后proc mean的clm选项能求统计量对应的置信上下限(注意也要取指数)
三:R语言的实现
一开始chatgpt给我的代码如下:
但是R算出的结果和SAS不一致
回到程序,注意红框中的qnorm,解释看的我有些不明白,我们查看SAS的proc means计算置信区间背后的统计公式,应用了t分布。所以问题应该就出在这个qnorm。
知道了SAS的proc means计算CI采用的是t分布,那就简单了,class数据集sex="F"有9个人,那么n-1就是8,我们去查t分布界值表(双侧)
我们将2.306带到R里面的程序中,z_score <- 2.306
bingo,看到算出来的置信区间跟SAS一样(分组为F),但是SEX=M的不一样,因为这组的观测是n=10,所以界值2.306并不适用于这组。所以平时了解点统计知识还是有用的。
其实R程序里面的那些代码,就是根据SASHELP里面的公式来的,现在知道问题在哪了,但是怎么解决我还不清楚,有时间研究一下这个norm,可能需要换成其他函数。
TA还问我是在哪学到这些东西
1)通过我的文章,可以看到我的思路,需要了解什么,首选SASHELP,最权威最详细,但是太枯燥难懂,而且有些选项不是都有对应的例子
2)第二就是sas community,平时有什么问题,搜索关键字,看下别人是不是有相同的问题,这个也解决了我不少问题
3)现在类似chatgpt的应用和软件太多了,有什么问题,直接把问题抛给它们,大多数都能给出正确的答案,但是不能太依赖,有些给的看似是正确答案,但是却是错的或者有误,还是需要自己认真推敲。