~ubuntu-branches/ubuntu/trusty/mricron/trusty

procedure FirthAnalyzeNoThread(lnCond, lnCrit,lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImagesCount : integer; lPlankImg: bytep;lOutImgMn,lSymptomRA: SingleP;lOutImg: SingleRAp);

implementation

uses npmform, dialogs;

procedure VisualProg(lPos: Integer);

begin

MainForm.ProgressBar1.Position := lPos;

MainForm.Refresh;

end;

var

finalloglik: SingleP0;

KxKA1,KxKB1,KxKA,KxKB :TMatrix;

Kvec,Kvec1 : TVector;

Kveci,kVeci1 : TVectori;

betak,xbeta,y,pi,ustar,

XXx,XXXW2,XXFisher,XXcovs,XXXWPrime,

deltahalfs,deltat,delta,covs,x,Fisher,XW2,W,XWprime,Hprime,H,ustarmat,negx: TMatrix;

lBarX: TProgressBar;

lnCondx,lnCritx,lBarPosX,lnPermuteX,lThreadx,lThreadStartx,lThreadEndx,lStartVoxx,lVoxPerPlankx,lImagesCountx : integer;

lPlankImgx: byteP;lOutImgMnx,lSymptomRAx: SingleP;

lOutImgX: SingleRAp;

procedure logistfx (xin: SingleP; var lZvals: SingleP0; numSubj,numCond: integer; lComputeIntercept: boolean);

//todo zero output incase exit

//yin = 1..numSubj binary 0/1 values

//xin = numSubj*numCond predictors

//Chivals = 0..numCond p-values - the 0th Khi-value is the intercept

// [0th value will not be computed if ; lComputeIntercept= false]

label

123,666;

const

maxit = 25;

maxhs = 5;

epsilon = 0.0001;

maxstep = 10;

var

SumY0,SumY1,mx, beta0,loglik,loglikold: double;

sumy, n, i,j, k, iter,halfs,lCond,dropCond: integer;

variability,firth: boolean;

procedure crossprodustar;

var

inc,row: integer;

begin

for row := 1 to k do begin

ustarmat[row,1] := 0;

for inc := 1 to ustar.r do

ustarmat[row,1] := ustarmat[row,1] + (x[row,inc]*ustar[inc,1]);

end;

procedure Diag2Vec;

var

inc: integer;

begin

for inc := 1 to pi.r do

ustar[inc,1] := ustar[inc,1]+ H[inc,inc]*(0.5-pi[inc,1]);

end; //nested DiagP2

procedure DiagP2 (var W, P: TMatrix);

var

inc: integer;

begin

W.Zero;

for inc := 1 to P.r do

W[inc,inc] := Power((P[inc,1] * (1-P[inc,1])),0.5) ;

end; //nested DiagP2

procedure ComputeFisher;

begin

DiagP2(W,pi);

XW2.mult(x,W);

//XWPrime.copy( XW2);

//XWPrime.transpose;

XWPrime.transpose(XW2);

Fisher.mult(XW2,XWPrime);

covs.copy( Fisher);

covs.Invert2(KxKA,KxKB,Kvec,Kveci)

end; //nested computeFisher

procedure computedropdelta;

var

jinc,iinc,ii,jj: integer;

begin

DiagP2(W,pi);

XXXW2.mult(XXx,W);

//XXXWPrime.copy( XXXW2);

//XXXWPrime.transpose;

XXXWPrime.transpose(XXXW2);

XXFisher.mult(XXXW2,XXXWPrime);

XXcovs.copy( XXFisher);

100

//XXcovs.Invert;

101

XXcovs.Invert2(KxKA1,KxKB1,Kvec1,Kveci1);

102

covs.Zero;

103

ii := 0;

104

for iinc := 1 to (k) do begin

105

if iinc <> (dropCond+1) then begin //leave the specified column zeros...

106

inc(ii);

107

jj := 0;

108

for jinc := 1 to (k) do begin

109

if jinc <> (dropCond+1) then begin

110

inc(jj);

111

covs[iinc,jinc] := xxCovs[ii,jj];

112

end;

113

end;

114

end;

115

end;

116

end;

117

function firthpenalty: double;

118

begin

119

ComputeFisher;

120

//result := 0.5 * ln(abs(Fisher.det));

121

result := 0.5 * ln(abs(Fisher.Det2(KxKA,kVeci,kVec)));

122

end; //nested firthpenalty

123

function ComputeLogLik: double;

124

var

125

inc: integer;

126

lDenom: double;

127

begin

128

xbeta.mult(betak,negx);

129

for inc := 1 to n do begin

130

lDenom := (1 + exp( xbeta[inc,1]));

131

if lDenom = 0 then

132

showmessage('yikes')

133

else

134

pi[inc,1] := 1/lDenom;

135

end;

136

result := 0;

137

for inc := 1 to n do

138

if y[inc,1] = 1 then

139

//if pi[inc,1] <> 1 then

140

result := result+ln(pi[inc,1]);

141

for inc := 1 to n do

142

if y[inc,1] = 0 then

143

//if pi[inc,1] <> 1 then

144

result := result+ln(1-pi[inc,1]);

145

if firth then

146

result := result + firthpenalty;

147

end;//nested ComputeLogLik

148

begin

149

for i := 0 to (numCond) do

150

lZVals^[i] := 0; //

151

if (numSubj < 2) or (numCond < 1) then

152

exit;

153

//ensure there is some variability in the input data...

154

variability := false;

155

i := 1;

156

repeat

157

inc(i);

158

if xin^[i] <> xin^[1] then

159

variability := true;

160

until (i= (numSubj*numCond)) or (variability);

161

if not variability then

162

exit; //no variance in the regressors...

163

variability := false;

164

i := 1;

165

repeat

166

inc(i);

167

if y[i,1] <> y[1,1] then

168

variability := true;

169

until (i= (numSubj)) or (variability);

170

if not variability then

171

exit; //no variance in the dependent variable...

172

dropCond := -1; //initially compute full model, then compute effect of removing individual conditions

173

firth := true;

174

n := numSubj;

175

k := numCond + 1;

176

//get memory

177

//beta := TMatrix.Create(n,1);

178

//design our model

179

//first row = 1: ell samples have equal weight

180

for i := 1 to n do

181

x.M[1,i] := 1;

182

//next load model into x

183

iter := 0;

184

for j := 2 to k do

185

for i := 1 to n do begin

186

inc(iter);

187

x.M[j,i] := xin^[iter];

188

end;

189

//WriteMatrix('Observations',y);

190

//WriteMatrix('Model',x);

191

//negx is just sing-swapped - we will generate this as we use it a lot...

192

for j := 1 to k do

193

for i := 1 to n do begin

194

negx.M[j,i] := -x.M[j,i];

195

end;

196

//now start computations

197

sumy := 0;

198

for i := 1 to n do

199

sumy := sumy + round(y[i,1]);

200

if (sumy <= 0) or (sumy >= n) then begin

201

//serious error: no variability. This should have been detected earlier in the procedure when yin was tested for variability

202

goto 666;

203

end;

204

beta0 := ln((sumy/n)/(1 - sumy/n));//initial estimate

205

123: //go here for each dropcond

206

if DropCond >= 0 then begin

207

betak.Ones;

208

betak.mult( 0) //start with a null model... does not really make sense

209

end else begin

210

betak.zero;

211

betak[1,1] := (beta0);

212

end;

213

iter := 0;

214

if DropCond >= 0 then begin //drop one of the factors...

215

if dropCond <> 0 then begin//include intercept

216

for i := 1 to n do

217

XXx.M[1,i] := 1;

218

lCond := 1;

219

end else

220

lCond := 0;

221

for j := 1 to NumCond do begin

222

if j <> DropCond then begin

223

inc(lCond);

224

for i := 1 to n do

225

XXx.M[lCond,i] := x.M[j+1,i];

226

end; //if j <> dropCond

227

end;

228

end;//if lDropCond >= 0

229

loglik := ComputeLogLik;

230

repeat

231

inc(iter);

232

ComputeFisher;

233

HPrime.mult(XWPrime,covs);

234

H.mult(HPrime,XW2);

235

//WriteMatrix(covs);

236

ustar.Sub(y,pi);

237

if firth then

238

Diag2Vec;

239

crossprodustar;

240

if dropCond >= 0 then // model with dropped factor

241

computedropdelta;

242

deltat.mult(covs,ustarmat);

243

delta.transpose(deltat);

244

mx := delta.MatAbsMax/MaxStep;

245

if mx > 1 then

246

delta.mult(mx);//scale delta

247

betak.add(delta);

248

loglikold := loglik;

249

halfs := 1;

250

while halfs <= maxhs do begin // Half-Steps

251

//fx(iter,halfs,loglik);

252

loglik := ComputeLogLik;

253

deltahalfs.mult(delta,power(2,-halfs));

254

betak.sub(deltahalfs);

255

if (loglik > loglikold) then

256

break;

257

inc(halfs);

258

end;

259

if delta.MatAbsMax <= epsilon then break;

260

until (iter >= maxit);

261

//fx(DropCond,loglik);

262

//done with this model - record model fit

263

if DropCond < 0 then

264

finalloglik^[k] := loglik //full model

265

else begin

266

finalloglik^[DropCond] := loglik; //model with a factor removed

267

end;

268

if DropCond < numCond then begin

269

inc(DropCond);

270

if (DropCond = 0) and (not lComputeIntercept) then //only compute intercept model if requested

271

inc(DropCond);

272

goto 123;

273

274

end;

275

//finally - results

276

277

//ResultsForm.Memo1.lines.add (inttostr(j)+' cases have Y=0, '+inttostr(n-j)+' cases have Y=1');

278

if lComputeIntercept then

279

J := 0

280

else

281

J := 1;

282

for i := J to (k-1) do begin

283

lZVals^[i] := abs(2*(finalloglik^[i]-finalloglik^[k]));

284

//find direction of effect - does a larger value of the IV predict more zeros or ones

285

lZVals^[i] := pNormalInv(ChiSq(lZVals^[i],1));

286

//we have now computed a Z scores - but Chi is one tailed, so all Z > 0... lets check direction

287

Sumy0 := 0;

288

Sumy1 := 0;

289

for iter := 1 to n do begin

290

if y[iter,1] = 0 then

291

Sumy0 := Sumy0 + x.M[i+1,iter] //+1: M indexed from 1, ZVal indexed from 0

292

else

293

Sumy1 := Sumy1 + x.M[i+1,iter]; //+1 M indexed from 1

294

end;

295

//compute means

296

Sumy1 := Sumy1/sumy;

297

Sumy0 := Sumy0/(n-sumy);

298

if Sumy0 < Sumy1 then //negative z-scores: damage here predicts performance is BETTER

299

lZVals^[i] := -lZVals^[i];

300

end;

301

(*if lComputeIntercept then //intercept is the 0th value

302

lChiVals[0] := abs(2*(finalloglik[0]-finalloglik[k]));

303

for i := 1 to (k-1) do //k-1 as this is indexed from 0

304

lChiVals[i] := abs(2*(finalloglik[i]-finalloglik[k])); *)

305

306

666:

307

end;

308

309

//FirthAnalyzeNoThread (lnCond,lnCrit, lnPermute,1,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,lPlankImg,lOutImgSum,lSymptomRA,lOutImg);

310

311

312

//calls logistf (yin,xin: SingleP; var lChivals: SingleP0; numSubj,numCond: integer);

313

label

314

666;

315

const

316

knPrevPattern = 10;

317

var

318

319

320

lPrevPatternRA: array[1..knPrevPattern] of TLesionPattern;

321

lPattern: TLesionPattern;

322

lObs: Bytep;

323

lPrevZVals: array [1..knPrevPattern] of SingleP0;

324

lZVals: SingleP0;

325

lPatternPos,lC,lnLesion,lPosPct,lPos,lPos2,lPos2Offset,lnCritLocal,n,k: integer;

326

begin //statthread

327

328

lnCritLocal := lnCrit;

329

if lnCritLocal < 1 then

330

lnCritLocal := 1;

331

Getmem(lObs,lImagesCount*sizeof(byte));

332

Getmem(lZVals,(lnCond+1)*sizeof(single));

333

for lPos := 1 to knPrevPattern do

334

Getmem(lPrevZVals[lPos],(lnCond+1)*sizeof(single));

335

n := lImagesCount;

336

k := lnCond + 1;

337

y := TMatrix.Create(n,1);

338

GetMem(finalloglik,(k+1)*sizeof(single));//finalloglik := TVector.Create(k+1);

339

x := TMatrix.Create (k, n);

340

betak:=TMatrix.Create(1,k);

341

covs:=TMatrix.Create(k,k);

342

delta:=TMatrix.Create(1,k);

343

deltahalfs:=TMatrix.Create(1,k);

344

deltat:=TMatrix.Create(k,1);

345

Fisher:=TMatrix.Create(k,k);

346

H:=TMatrix.Create(n,n);

347

HPrime:=TMatrix.Create(n,k);

348

negx:=TMatrix.Create(k,n);

349

pi:=TMatrix.Create(n,1);

350

ustar:=TMatrix.Create(n,1);

351

ustarmat:=TMatrix.create(k,1);

352

W:=TMatrix.Create(n,n);

353

xbeta:=TMatrix.Create(1,n);

354

XW2:=TMatrix.Create(k,n);

355

//XWPrime:=TMatrix.Create(k,n);

356

XWPrime:=TMatrix.Create(n,k);

357

XXcovs:=TMatrix.Create(k-1,k-1);

358

XXFisher:=TMatrix.Create(k-1,k-1);

359

XXx:=TMatrix.Create(k-1,n);

360

XXXW2:=TMatrix.Create(k-1,n);

361

//XXXWPrime:=TMatrix.Create(k-1,n);

362

XXXWPrime := TMatrix.Create ( n, k-1);

363

KxKA := TMatrix.Create(k,k);

364

KxKB := TMatrix.Create(k,k);

365

Kvec := TVector.Create(k);

366

Kveci := TVectori.Create(k);

367

KxKA1 := TMatrix.Create(k-1,k-1);

368

KxKB1 := TMatrix.Create(k-1,k-1);

369

Kvec1 := TVector.Create(k-1);

370

Kveci1 := TVectori.Create(k-1);

371

372

lPosPct := (lThreadEnd-lThreadStart) div 100;

373

for lPatternPos := 1 to knPrevPattern do

374

lPrevPatternRA[lPatternPos] := EmptyOrder;

375

lPatternPos := 1;

376

for lPos2 := lThreadStart to lThreadEnd do begin

377

if (lThread = 1) and ((lPos2 mod lPosPct) = 0) then

378

VisualProg(round((lPos2/(lThreadEnd-lThreadStart))*100));

379

lPos2Offset := lPos2+lStartVox-1;

380

lnLesion := 0;

381

382

for lPos := 1 to lImagesCount do begin

383

if lPlankImg^[((lPos-1)* lVoxPerPlank)+lPos2] = 0 then begin

384

//no lesion

385

y[lPos,1] := 0;

386

lObs^[lPos] := 0;

387

end else begin

388

//lesion

389

inc(lnLesion);

390

lObs^[lPos] := 1;

391

y[lPos,1] := 1; //note: lObs indexed from zero!

392

end;

393

end;

394

lOutImgMn^[lPos2Offset] := lnLesion;///lImages.Count;

395

if (lnLesion >= lnCritLocal) and (lnLesion < lImagesCount) then begin

396

lPattern := SetOrderX (lObs,lImagesCount);

397

lPos := 1;

398

while (lPos <= knPrevPattern) and not (SameOrder(lPattern,lPrevPatternRA[lPos],lImagesCount)) do

399

inc(lPos);

400

if SameOrder(lPattern,lPrevPatternRA[lPos],lImagesCount) then begin

401

inc(gnVoxTestedRA[lThread]);

402

//logistfx(lObs,lSymptomRA, lZvals, lImagesCount,lnCond,false);

403

for lC := 1 to lnCond do

404

lOutImg^[lC]^[lPos2Offset] := lPrevZvals[lPos]^[lC];

405

end else begin //new pattern - need to compute

406

inc(gnVoxTestedRA[lThread]);

407

logistfx(lSymptomRA, lZvals, lImagesCount,lnCond,false);

408

for lC := 1 to lnCond do

409

lOutImg^[lC]^[lPos2Offset] := lZvals^[lC];

410

lPrevPatternRA[lPatternPos] := lPattern;

411

for lC := 1 to lnCond do

412

lPrevZVals[lPatternPos]^[lC] := lZvals^[lC];

413

inc(lPatternPos);

414

if lPatternPos > knPrevPattern then

415

lPatternPos := 1;

416

417

end; //new pattern

418

end; //nlesion > nCritical

419

420

end; //for each voxel

421

//gMat := false;

422

666:

423

freemem(lObs);

424

for lPos := 1 to knPrevPattern do

425

freemem(lPrevZVals[lPos]);

426

freemem(lZVals);

427

428

y.free;

429

x.free;

430

betak.free;

431

covs.free;

432

delta.free;

433

deltahalfs.free;

434

deltat.free;

435

Fisher.free;

436

H.free;

437

HPrime.free;

438

negx.free;

439

pi.free;

440

ustar.free;

441

ustarmat.Free;

442

W.free;

443

xbeta.free;

444

XW2.free;

445

XWPrime.free;

446

XXcovs.free;

447

XXFisher.free;

448

XXx.free;

449

XXXW2.free;

450

XXXWPrime.free;

451

KxKA.free;

452

KxKB.free;

453

Kvec.free;

454

Kveci.free;

455

KxKA1.free;

456

KxKB1.free;

457

Kvec1.free;

458

Kveci1.free;

459

460

freemem(finalloglik);

461

462

end;

463

464

end.

465

b'\\ No newline at end of file'

Older »