R语言数据分析.zip

# PROCESS for R version 4.1.1 # Written by Andrew F. Hayes # www.afhayes.com # www.processmacro.org # Copyright 2012-2022 by Andrew F. Hayes ALL RIGHTS RESERVED # PROCESS workshop schedule at http://haskayne.ucalgary.ca/CCRAM # # Distribution of this code in any form except through processmacro.org # is prohibited without the permission of the copyright holder, as is # distribution after modification. # # THIS SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT # IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM # DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT # OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE # USE OF THIS SOFTWARE IMPLIES AGREEMENT WITH THESE TERMS # # To activate, run this script. It may take a few minutes. # Patience is a virtue. This will produce a function called process. # The other functions it creates cannot be accessed by the user # but are used by the master process function. process.bcboot3<-function(databcbt,estmte,xp2,badend,priorlo,priorhi) { databcbt<-as.matrix(sort(databcbt)) badlo<-0;badhi<-0 pv<-matrix(as.numeric(databcbt < estmte)); pv<-sum(pv)/nrow(databcbt);ppv<-pv; if (pv > .5){ppv<-(1-pv)} y5<-sqrt(-2*log(ppv)) p0<-(-.322232431088);p1<-(-1);p2<-(-.342242088547);p3<-(-.0204231210245) p4<-(-.0000453642210148);q0<-(.0993484626060);q1<-(.588581570495) q2<-(.531103462366);q3<-(.103537752850);q4<-(.0038560700634) xp<-y5+((((y5*p4+p3)*y5+p2)*y5+p1)*y5+p0)/((((y5*q4+q3)*y5+q2)*y5+q1)*y5+q0) if (pv <= .5){xp<-(-xp)} cilow<-round(nrow(databcbt)*pnorm(2*xp-xp2)) cihigh<-trunc(nrow(databcbt)*pnorm(2*xp+xp2))+1 if (cilow < 1){cilow<-1;booterr<-1;badlo<-1} if (cihigh > nrow(databcbt)){cihigh<-nrow(databcbt);booterr<-1;badhi<-1} llcit<-databcbt[cilow,1] ulcit<-databcbt[cihigh,1] if ((badlo==1) & (llcit != priorlo)){priorlo<-llcit;badend<-c(badend,llcit)} if ((badhi==1) & (ulcit != priorhi)){priorhi<-ulcit;badend<-c(badend,ulcit)} bootse<-sd(databcbt) cires<-as.matrix(c(bootse,llcit,ulcit)) cires<-list(cires,badend,priorlo,priorhi) return(cires) } process.pboot3<-function(databcbt,lcval,hcval) { databcbt<-as.matrix(sort(databcbt)) llcit<-databcbt[lcval,1] ulcit<-databcbt[hcval,1] bootse<-sd(databcbt) cires<-as.matrix(c(bootse,llcit,ulcit)) return(cires) } process.llrtest3<-function(lm,y,x,b,basemod,iterate,converge) { lm<-as.matrix(lm) btemphld<-b llrdat<-matrix(-999,nrow(x),(nrow(lm)-sum(lm))) llrdf<-ncol(x)-ncol(llrdat) llrcnt<-0 for (llri in (1:nrow(lm))) { if (lm[llri,1]==0){llrcnt<-llrcnt+1;llrdat[,llrcnt]<-x[,llri]} } LL2<-process.modelest(y,llrdat,2,0,xp2,5,iterate,converge) b<-btemphld pvchi<-(1-pchisq((LL2-basemod),df=llrdf)) fresult<-cbind((LL2-basemod),llrdf,pvchi) return(fresult) } process.describ3<-function(descdatf,type=0,quantle=1) { desctmp<-matrix(-999,(8-(4*type)),ncol(descdatf)) # mean, sd, min, max, 16th, 50th, 84th, dich toggle for (jd in c(1:ncol(descdatf))) { descdat<-descdatf[,jd] #get the mean, sd, minimum, and maximum */ desctmp[1,jd]<-mean(descdat) desctmp[2,jd]<-sd(descdat) desctmp[3,jd]<-min(descdat) desctmp[4,jd]<-max(descdat) if (type==0) { minwarn<-0;maxwarn<-0 tmp=as.numeric(descdat==desctmp[3,jd])+as.numeric(descdat==desctmp[4,jd]) desctmp[8,jd]<-as.numeric(sum(tmp)==length(tmp)) if (desctmp[3,jd]==desctmp[4,jd]){desctmp[8,jd]<-2} descdat<-matrix(sort(descdat)) decval<-c(.16,.5,.84) for (kd in c(1:3)) { low<-trunc(decval[kd]*(length(descdat)+1)) lowdec<-decval[kd]*(length(descdat)+1)-low value<-descdat[low,1]+(descdat[(low+1),1]-descdat[low,1])*lowdec desctmp[(4+kd),jd]<-value } mnotev<-(1) modvals<-matrix(desctmp[5:7,],ncol=ncol(descdatf)) if (quantle != 1) { desctmp[5,jd]<-desctmp[1,jd]-desctmp[2,jd] desctmp[6,jd]<-desctmp[1,jd] desctmp[7,jd]<-desctmp[1,jd]+desctmp[2,jd] modvals<-matrix(desctmp[5:7,],ncol=ncol(descdatf)) mnotev<-(2) if (modvals[1,1] < desctmp[3,1]){modvals[1,1]<-desctmp[3,1];minwarn<-1} if (modvals[3,1] > desctmp[4,1]){modvals[3,1]<-desctmp[4,1];maxwarn<-1} } if (desctmp[8,jd]==1) {modvals<-matrix(c(desctmp[3,1],desctmp[4,1])) mnotev<-0;minwarn<-0;maxwarn<-0 } descrtrn<-list(desctmp,modvals,minwarn,maxwarn,mnotev) } } if (type==1) {descrtrn<-list(desctmp)} return(descrtrn) } process.ftest3<-function(lm,bcoef,cv=0,chr=0,brsq=0,skip=0,y,x) { lmat2<-as.matrix(lm) y<-as.matrix(y) x<-as.matrix(x) n<-nrow(y) if (skip==0) { lmat2<-as.matrix(diag(as.numeric(lm))) lmat3<-matrix(0,nrow(lmat2),1) for (flp in c(1:ncol(lmat2))) { if (sum(lmat2[,flp])==1) {lmat3<-cbind(as.matrix(lmat3),as.matrix(lmat2[,flp]))} } lmat2<-as.matrix(lmat3[,2:ncol(lmat3)]) } fratio<-(t(t(lmat2)%*%bcoef)%*%solve(t(lmat2)%*%cv%*%lmat2)%*%((t(lmat2)%*%bcoef)))/ncol(lmat2) pfr<-(1-pf(fratio,ncol(lmat2),(n-nrow(bcoef)))) fresult<-matrix(c(fratio,ncol(lmat2),(n-nrow(bcoef)),pfr),ncol=4) if (chr==1) { lmat3<-as.matrix(1-rowSums(lmat2)) xfm<-matrix(0,n,sum(lmat3)) flpc<-1 for (flp in (1:nrow(lmat3))) { if (lmat3[flp,1]==1){xfm[,flpc]=x[,flp];flpc<-flpc+1} } bfm<-solve(t(xfm)%*%xfm)%*%t(xfm)%*%y resid<-y-(xfm%*%bfm) sstotal<-t(y-(sum(y)/n))%*%(y-(sum(y)/n)) ssresid<-t(resid)%*%resid rsqch<-as.numeric(brsq)-((sstotal-ssresid)/sstotal) fresult<-matrix(c(rsqch,fresult),ncol=5) } ftestout<-as.matrix(fresult) return(ftestout) } #type1 for ols, type2 for logistic LLR, type3 for logistic bootstrapping process.modelest<-function(y,x,type,full,xp2,hc=5,iterate=100,converge=.00001) { if (type==1) { invxtx<-solve(t(x)%*%x) b<-invxtx%*%t(x)%*%y modres<-b if (full==1) { n1<-nrow(x) dfres<-(n1-(ncol(x))) sstotal<-t(y-(sum(y)/n1))%*%(y-(sum(y)/n1)) resid=y-x%*%b ssresid<-sum(t(resid)%*%resid) r2<-(sstotal-ssresid)/sstotal adjr2<-(1-((1-r2)*(n1-1)/(dfres))) mse<-ssresid/(n1-ncol(x)) #HC covariance matrix varb<-mse*invxtx k3<-ncol(x) xhc<-0 if (hc != 5) { xhc<-x hat<-matrix(xhc[,1]) for (i3 in c(1:nrow(xhc))) { xhcm<-matrix(xhc[i3,]) hat[i3,1]<-t(xhcm)%*%invxtx%*%xhcm } if ((hc==0) | (hc==1)) { for (i3 in c(1:k3)){xhc[,i3]<-xhc[,i3]*resid} } if ((hc==3) | (hc==2)) { for (i3 in c(1:k3)) {xhc[,i3]<-(resid/(1-hat)^(1/(4-hc)))*xhc[,i3]} } if (hc==4) { hcmn<-matrix(4,n1,2);hcmn[,2]<-(n1*hat)/k3 minr<-apply(hcmn,1,FUN=min) for (i3 in c(1:k3)) { xhc[,i3]<-(resid/(1-hat)^(minr/2))*xhc[,i3] } } varb<-(invxtx%*%t(xhc)%*%xhc%*%invxtx) if (hc==1){varb<-(n1/(n1-ncol(x)))*varb} } seb<-sqrt(diag(varb)) trat<-b/seb p<-2*pt(-abs(trat),df=dfres) tval<-sqrt(dfres* (exp((dfres-(5/6))*((xp2/(dfres-(2/3)+(.11/dfres)))*(xp2/(dfres-(2/3)+(.11/dfres)))))-1)) modres<-matrix(c(modres,seb, trat,p,(b-tval*seb),(b+tval*seb)),ncol=6) modresl<-t(matrix(c("coeff","hclab","t","p","LLCI","ULCI"))) lmat<-diag(ncol(x));lmat<-lmat[,2:ncol(lmat)] fratio<-(t(t(lmat)%*%b)%*%solve(t(lmat)%*%varb%*%lmat)%*%((t(lmat)%*%b)))/(ncol(x)-1) pfr<-1-pf(fratio,(ncol(x)-1),dfres) modsum=matrix(c(sqrt(r2),r2,mse,fratio,(ncol(x)-1),dfres,pfr)) modsuml=matrix(c("R","R-sq","MSE","hcflab","df1","df2", "p")) modretrn<-list(modres,modresl,modsum,modsuml,b,varb,tval,resid) return(modretrn) } if (full==0){return(modres)} } #for logistic Y model if ((type==2) | (type==3)) { xlp<-x;ylp<-as.matrix(y) pt2<-matrix((sum(ylp)/nrow(ylp)),nrow(ylp),1) if ((type==2)|(type==3)) {LL3<-(ylp*log(pt2)+(1-ylp)*log(1-pt2))} LL3<-(-2*sum(LL3)) bt1<-matrix(0,ncol(xlp),1);LL1<-0 pt1<-matrix(0.5,nrow(ylp),1);pt1lp<-pt1 for (jjj in (1:iterate)) { xlptmp<-t(xlp) vecprb<-(pt1lp*(1-