Nested pure functions can make this work. I write the integrand as a nested pure function, with the outer function taking, as an argument, the variable of integration, and the inner pure function taking the solution of the underlying equation and using it to produce an explicit list of the resulting objects of interest.

Clear[x, equation, startVals, modelOutcomes]; startVals = Range[1/10, 9/10, 1/10]; equation[x_?NumericQ, param_] := Beta[2, 2] - Beta[x, 2, 2] - x*(1 - x)*param; xSoln[param_?NumericQ] := Min[Part[FindRoot[equation[xVal, param], {xVal, #, 0, 1}] & /@ startVals, All, 1, 2]]; result1 = Beta[#, 2, 2]/Beta[2, 2] &; result2 = Beta[#, 3, 2]/Beta[#, 2, 2] &; modelOutcomes = {result1[#], result2[#]} &[xSoln[#]] &; mu = 1/10; sigma = 1/1000; NExpectation[modelOutcomes[z1], z1 \[Distributed] LogNormalDistribution[mu, sigma]] (*{0.0374898, 0.0550367}*) 

I would prefer not to use nested pure functions, because I find them more difficult to read, but I don't know how else to accomplish this.

We can write the integrand as a nested pure function, with the outer function taking, as an argument, the variable of integration, and the inner pure function taking the solution of the underlying equation and using it to produce an explicit list of the resulting objects of interest.

Clear[x, equation, startVals, modelOutcomes]; startVals = Range[1/10, 9/10, 1/10]; equation[x_?NumericQ, param_] := Beta[2, 2] - Beta[x, 2, 2] - x*(1 - x)*param; xSoln[param_?NumericQ] := Min[Part[FindRoot[equation[xVal, param], {xVal, #, 0, 1}] & /@ startVals, All, 1, 2]]; result1 = Beta[#, 2, 2]/Beta[2, 2] &; result2 = Beta[#, 3, 2]/Beta[#, 2, 2] &; modelOutcomes = {result1[#], result2[#]} &[xSoln[#]] &; mu = 1/10; sigma = 1/1000; NExpectation[modelOutcomes[z1], z1 \[Distributed] LogNormalDistribution[mu, sigma]] (*{0.0374898, 0.0550367}*) 

This approach has a disadvantage, which is that NIntegrate is essentially separately integrating each element of the integrand vector. I can't demonstrate this with NExpectation, because it doesn't have the EvaluationMonitor option, but I can demonstrate it with NIntegrate. So for demonstration purposes I'll use a different function, but using the same nested pure functions approach to the vector integrand

Clear[g, f, samplePointsList]; g[y_] := y^2; f := {Sin[#], Cos[#]} &[g[#]] & 

Borrowing an approach from this question, I collect the evaluation points used by NIntegrate, and plot the point on the horizontal axis with its position in the order of evaluated points on the vertical axis,

samplePointsList = Reap[NIntegrate[f[x], {x, 0, 1}, EvaluationMonitor :> Sow[x] ]][[2, 1]]; ListPlot[Transpose[{samplePointsList, Range[Length[samplePointsList]]}]] 

It looks like the evaluation happens twice at essentially the same set of points. What fraction of the evaluation points are unique?

N@Length[Union[samplePointsList]]/Length[samplePointsList] (*0.505792*) 

So this approach facilitates the use of NIntegrate / NExpectation, but doesn't deal with the duplicated evaluation problem.