title: Decorrelating two image maps

text: I am stuck on a research problem and do not know what kind of methods to use. I hope people in this forum can give me some good ideas, I've always had good brainstorming sessions here. My problem is that I have a laser and a stack of thin films made up of three different materials. If I shoot the laser at these thin films I get a circular image of the reflected light from the laser scattering on the films. However, the fourth layer is added and I want to measure the change in thickness of the fourth layer, similarly by shooting a laser at this stack. The other three bottom layers are "nuisance" parameters, but their thickness variations are also captured in the pupil image. How do I filter out the observed pupil image so that the information of the three base stacks are eliminated and only the variations of the fourth critical stack are captured? I have 1000 sets (1000 for 3-stack and 1000 for 4-stack) such that for each data point, the 3-stack and 4-stack films have gone through identical process conditions and there is an image for both these stacks. I've thought of maybe signal transforms? ------- edit -------- Each image pixel of the pupil of the film stack is a nonlinear function of the thicknesses but the gradients at each parameter (as well as the image at a certain parameter value) point can be calculated via an optical simulator. $pix_i(i = 1...n) = f(h_1,h_2,h_3) + \epsilon_i$ where $pix_i$ is a pixel point of the acquired image and $h1,h2,h3$ are the film thicknesses of a 3-film stack. This can now be linearized around $h10,h20,h30$ to $pix_i(i=1...n) \sim f(h_{10},h_{20},h_{30}) + \frac{df}{dh_i}(h_i = h_{i0}) + \epsilon_i$


title: Relations between probabilities of "almost" independent random variables

text: Let $X$ and $Y$ be two random variables, such that the (average) mutual information is very small: $$ 0 \le I(X;Y) \le \epsilon \ll 1$$ In this case, we say that $X$ and $Y$ are _almost_ independent. Now, can we deduce something like: $$\forall x,y \quad \Pr[X=x,Y=y] \le \Pr[X=x]\cdot\Pr[Y=y] + \delta$$ with $\delta$ being a (possibly negative) function of $\epsilon$? * * * If the general case above cannot be proven, maybe we can use the following special case: Let $X$ and $Y$ be **discrete** random variables, taking values from a **finite** set $D$ with $|D|=2^n$. Moreover, let $\epsilon$ be exponentially small in $n$, say $\epsilon = 2^{-n/2}$. Can we now obtain $\delta$ as a function of $\epsilon$?


title: How should we convert sports results data to perform a valid logistical regression?

text: Say we want to perform a logistical regression analysis (although my question pertains to regressions in general) on sports results to determine the effects of various factors on who wins and who loses. We have the background information we want on the teams and players and now just need a random sample. So we decide to take the published results over the past couple of years as our sample. The sample we collect is in the following form: Result, Team 1, Team 2, ... The result is always 0-1 or 1-0 (no draws). We can start preparing the data by converting Result into a binary variable: Result = 1 if Team 1 wins, = 0 if Team 2 wins. The problem is that this doesn't give us a valid regression. The reason will take a bit of explaining. Say one of our observations is: Result = 1; Team 1 = Man.U.; Team 2 = Chelsea This observation can be rewritten: Result = 0; Team 1 = Chelsea; Team 2 = Man.U. And it is the exact same observation and all the information is still the same and perfectly correct. And this actually changes the results of our regression! One quick way to prove this is to consider what happens if we rewrite all of the observations so that Team 1 always wins. Then our dependent variable will always be Result = 1. Thus Var(Result) = 0 and the estimates for our parameters will all be 0 (except for the constant, of course). If we flip half of the observations so that half the time Result = 1 and half the time Result = 0 and we run the regression on that, we will get non-zero estimates for our parameters. This bothers me because we are regressing the same data but getting wildly different results based on the order the teams are written in. If our results can change based on the order we decided to put the teams down when recording our observations, then they can't be valid. So what is the best way to prepare this data for analysis so that we can get valid results?


title: When does quantile regression produce biased coefficients (if ever)?

text: It is easy to show using matrix algebra when least squares will produce bias. \begin{equation} \begin{split} \text{E}[B]& = \text{E}[(X'X)^{-1}]\times\text{E}[X'Y] \\ & = \text{E}[(X'X)^{-1}]\times\text{E}[X'(XB +\epsilon)] \\ & = \text{E}[(X'X)^{-1}]\times\text{E}[X'XB] + \text{E}[(X'X)^{-1}]\times\text{E}[X'\epsilon] \\ & = B + \text{E}[(X'X)^{-1}]\times\text{E}[X'\epsilon] \\ \end{split} \end{equation} Bias is introduced in the last term when there is correlation between X and e. My question is, how do we know when $\text{E}[B]$ for the LAD estimator will be biased? Will bias arise when there is correlation between $X$ and $\epsilon$, as in least squares? Or is quantile regression robust to correlation between $X$ and $\epsilon$? I'm guessing it can't be demonstrated using matrix algebra because the LAD estimator is the following: \begin{align} B(\tau) = \operatorname{argmin} \text{E}[\rho(Y_i - X'B)] \end{align} and the LAD estimator is not calculated using linear algebra (?). If that's right, then how can we demonstrate when quantile regression will produce bias? Use a monte carlo simulation?


title: Does "Inference" include estimation or only testing?

text: Does the term "statistical inference" include only hypothesis testing or does it also include point estimation, interval estimation etc. Authoritative references will be greatly appreciated.


title: Using MatchIt to match groups in a retrospective analysis

text: I am interested in using the R package MatchIt to preprocess my data as to obtain matched groups based on a predefined treatment variable. However I am facing a few issues. The first issue is that my data contains many covariates (approximately 80) in addition to the treatment variable. Many of these covariates have missing values, and it seems that MatchIt does not support missing variables. For those of you who have used MatchIt, what is your approach to solving this problem? Do you use imputation methods (e.g. the package `mice`) avoid NAs? Or do you perform variable selection to choose the most significant variables (e.g. lasso using `glmnet`) to do dimension reduction +/- imputation? The second issue is that I have many observations (n = 23000) with about 200 in my treatment group and the remaining 22800 in my control group. In this situation which method in MatchIt would be best suited (also taking into consideration the other aspects of my data as I have described above)? As always, many thanks for the help.


title: How to get distribution of sum of dependent bernoulli variables

text: I have $N$ Bernoulli variables, $X_1,...,X_N$ and $X_i\sim B(1, \pi_i)$, $\pi$ is known for each $X_i$, and $Y=X_1+...+X_N$, now I need to get the destribution of $Y$. If $X_i$ and $X_j$ are independent when $i\ne j$, then I can use the simulation:               1. Generate X1, ..., XN via their distribution, and then get the value of Y;     2. Repeat step 1 for 10000 times, and then I can get Y1, ..., Y10000, so I can know the distribution of Y.      But now $X_i$ and $X_j$ are dependent, so I also need to take into account the correlation, assuming that $\text{cor}(X_i, X_j)=0.2$ when $i\ne j$, how can I insert the correlation to the simulation? Or get the distribution of $Y$ via other ways?


title: Combining evidence from continuous classification

text: I am trying to predict expert ratings. They were asked to assign confidence among five possibilities: variables `t1` through `t5`. I have a training set for which text allows me to predict these probabilities. Is there something that I'm not considering for a model of the kind: `1 = t1 + t2 + t3 + t4 + t5`? If I find several classification methods, how could I combine the evidence? The averages of the training set t-variables: `0.2, 0.05, 0.3, 0.1, 0.35` in order. Let's assume I predict using two methods for a single case:               method 1: 0.1, 0.1, 0.4, 0.1, 0.3     method 2: 0.15, 0.2, 0.15, 0.2, 0.3      I considered simply taking the mean, but that does not seem right. I also considered updating like: `1 - (1 - m1) * (1 - m2)`, but the means have to be taken into account, as evidence below the mean (e.g. combining 0.1 and 0.05 will go above 0.1, given that the mean is 0.2), while this should go in the opposite direction of the mean. What am I missing here? Any feedback is welcome. **EDIT:** Another thing I considered is to calculate effect sizes and do something with that, perhaps shrink the lower ones.


title: Cross validation and ordinal logistic regression

text: I am trying to understand cross-validation for ordinal logistic regression. The aim of the game is to validate the model used in an analysis... I first construct a toy data set:               set.seed(1)     N <- 10000     # predictors     x1 <- runif(N)     x2 <- runif(N)     x3 <- runif(N)          # coeffs in the model     a <- c(-2,-1)     x <- -x1+2*x2+x3          # P( y ≤ i ) is given by logit^{-1} ( a[i]+x )     p <- outer(a,x, function(a,x) 1/(1+exp(-a-x)) )          # computing the probabilities of each category     q <- 1 - p[2,]     p[2,] <- p[2,] - p[1,];     p <- rbind(p,q);          # outcome     y <- ordered( apply( p, 2, function(p) which(rmultinom(1,1,p)>0) ) )       Now, I fit the model it using `lrm` in the package `rms`.               require("rms")     fit <- lrm(y~x1+x2+x3, x=TRUE,y=TRUE)               > fit          Logistic Regression Model          lrm(formula = y ~ x1 + x2 + x3, x = TRUE, y = TRUE)                                Model Likelihood     Discrimination    Rank Discrim.                              Ratio Test            Indexes          Indexes        Obs         10000    LR chi2    1165.46    R2       0.126    C       0.664      1           2837    d.f.             3    g        0.779    Dxy     0.328      2           2126    Pr(> chi2) <0.0001    gr       2.178    gamma   0.329      3           5037                          gp       0.147    tau-a   0.203     max |deriv| 4e-10                          Brier    0.187                                Coef    S.E.   Wald Z Pr(>|Z|)     y>=2  2.1048 0.0656  32.06 <0.0001     y>=3  1.0997 0.0630  17.45 <0.0001     x1    0.8157 0.0675  12.09 <0.0001     x2   -1.9790 0.0701 -28.21 <0.0001     x3   -1.0095 0.0687 -14.68 <0.0001      I understand the second part of the outcome: the coefficients I put in the model are here (it is almost perfect with `N = 100000`). The sign is reversed because in my model I used the coeffs to compute the odds of being $\le 1$, and $\le 2$, here it’s the other way, I think there’s not much issues there. However **I don’t understand the discrimination and rank discrimination indexes**. Can you help me?! Some pointers? Things are worse when we turn to cross validation...               > validate(fit, method="cross")               index.orig training    test optimism index.corrected  n     Dxy           0.3278   0.3278  0.3290  -0.0012          0.3291 40     R2            0.1260   0.1260  0.1313  -0.0053          0.1313 40     Intercept     0.0000   0.0000 -0.0072   0.0072         -0.0072 40     Slope         1.0000   1.0000  1.0201  -0.0201          1.0201 40     Emax          0.0000   0.0000  0.0056   0.0056          0.0056 40     D             0.1164   0.1165  0.1186  -0.0021          0.1186 40     U            -0.0002  -0.0002 -0.8323   0.8321         -0.8323 40     Q             0.1166   0.1167  0.9509  -0.8342          0.9509 40     B             0.1865   0.1865  0.1867  -0.0001          0.1867 40     g             0.7786   0.7786  0.7928  -0.0142          0.7928 40     gp            0.1472   0.1472  0.1478  -0.0007          0.1478 40      Mmffff? **What’s this? How do I interpret this?** The man page gives few explaination, I don’t have access to this paper... and I feel overwhelmed by an ocean of complexity. Please help!


title: Understanding Bayesian Predictive Distributions

text: I'm taking an Intro to Bayes course and I'm having some difficulty understanding predictive distributions. I understand why they are useful and I'm familiar with the definition, but there are some things I don't quite understand. **1) How to get the right predictive distribution for a vector of new observations** Suppose that we have built a sampling model $p(y_i | \theta)$ for the data and a prior $p(\theta)$. Assume that the observations $y_i$ are conditionally independent given $\theta$. We have observed some data $\mathcal{D} = \{y_1, y_2, \, ... \, , y_k\}$, and we update our prior $p(\theta)$ to the posterior $p(\theta | \mathcal{D})$. If we wanted to predict a vector of new observations $\mathcal{N} = \{\tilde{y}_1, \tilde{y}_2, \, ... \, , \tilde{y}_n\}$, I think we should try to get the posterior predictive using this formula $$ p(\mathcal{N} | \mathcal{D}) = \int p(\theta | \mathcal{D}) p ( \mathcal{N} | \theta) \, \mathrm{d} \theta = \int p(\theta | \mathcal{D}) \prod_{i=1}^n p(\tilde{y}_i | \theta) \, \mathrm{d} \theta, $$ which is not equal to $$ \prod_{i=1}^n \int p(\theta | \mathcal{D}) p(\tilde{y}_i | \theta) \, \mathrm{d} \theta, $$ so the predicted observations are not independent, right? Say that $\theta | \mathcal{D} \sim$ Beta($a,b$) and $p(y_i | \theta) \sim$ Binomial($n, \theta$) for a fixed $n$. In this case, if I wanted to simulate 6 new $\tilde{y}$, if I understand this correctly, it would be wrong to simulate 6 draws independently from the Beta-Binomial distribution that corresponds to the posterior predictive for a single observation. Is this correct? I don't know how to interpret that the observations are not independent marginally, and I'm not sure I understand this correctly. **Simulating from posterior predictives** Many times when we simulate data from the posterior predictive we follow this scheme: For $b$ from 1 to $B$: 1) Sample $\theta^{(b)}$ from $p(\theta | \mathcal{D})$. 2) Then simulate new data $\mathcal{N}^{(b)}$ from $p(\mathcal{N} | \theta^{(b)})$. I don't quite know how to prove this scheme works, although it looks intuitive. Also, does this have a name? I tried to look up a justification and I tried different names, but I had no luck. **Thanks!**


title: ggplot2: geom_polygon with no fill

text: i hope, you dont need data for this problem, because i believe i am just making a stupid syntax error. The following code:               ggplot()+       geom_point(data=sites, aes(x=NMDS1, y=NMDS2, shape=group), colour="grey") +         geom_point(data=species, aes(x=NMDS1, y=NMDS2, color=phyla), size=3, shape=20) + scale_colour_manual(values=Pal1) +       geom_segment(data = BiPlotscores, aes(x = 0, xend = NMDS1, y= 0, yend = NMDS2),                    arrow = arrow(length = unit(0.25, "cm")), colour = "black") +       geom_text(data = BiPlotscores, aes(x = 1.1*NMDS1, y = 1.1*NMDS2, label = Parameters), size = 3) + coord_fixed()+       theme(panel.background = element_blank()) +       geom_polygon(data = hulls, aes(x=NMDS1, y=NMDS2, colour=phyla, alpha = 0.2))      leads to the following result: ![enter image description here](http://i.stack.imgur.com/wLv53.jpg) (This is not the final product :)). I would like to have the polygons unfilled, or very just neatly filled. I do not want them to be greyish, for sure. Fill doesnt do anything, and apparently fiddling with alpha doesnt change anything, either. Any ideas are superwelcome. Thank you very much! "Hulls" is coming from the following code (as found here somewhere):               #find hulls     library(plyr)     find_hull <- function(df) df[chull(df$NMDS1, df$NMDS2), ]     hulls <- ddply(species , "phyla", find_hull)


title: Intersubject concordance for ratings vs rankings

text: In my system I ask $n$ people to rate [1,5] $m$ objects (e.g. films). I'd like to get a measure of the inter-subject concordance, i.e. how much the users agreed overall, not just on a particular object but on all of the objects. I know that Kendall's W is a well-established measure of concordance on rankings, but in my case people are not ordering objects, but assigning a 1-5 score. From Kendall's perspective, this results in a lot of ties, making Kendall's W problematic. What is a suitable measure of inter-subject concordance for _ratings_ rather than for _rankings_? Thanks, Mulone


title: Convex vs. concave curve - measurement or index

text: I am interested in a measurement or index that will tell me if a curve is more convex or concave versus a straight line. This curve is sampled in n points and I know the x and y coordinates of these points. The straight line I am comparing this curve with is a line that goes through the first and last point of the curve. Let's suppose that the straight line is oriented bottom left to top right with a certain angle and if the curve is on the left side of the straight line or north of it (top) then it is convex and if it is on the right side or south (bottom) it is concave. If this curve never intersects the straight line – then I can imagine a way to maybe deal with this problem, but if it intersects the straight line then... I am at a loss. Any suggestions will be very much appreciated.


title: How do I improve my neural network stability?

text: I'm using the neuralnet in R to build a NN with 14 inputs and one output. I build/train the network several times using the same input training data and the same network architecture/settings. After each network is produced I use it on a stand alone set of test data to calculate some predicted values. I'm finding there is a large variance in each iteration of the predicted data, despite all the inputs (both the training data and test data) remaining the same each time I build the network. I understand that there will be differences in the weightings produced within the NN each time and that no two neural networks will be indentical, but what can I try to produce networks that are more consistent across each train, given the indentical data? Thanks in advance for any comments.


title: How to perform post-hoc comparison on interaction term with mixed-effects model?

text: I’m working on a data set in order to evaluate the impact of drying on sediment microbial activities. The objective is to determine if the impact of drying varies among sediment types and/or depth within the sediment. The experimental design is as follows:   * The first factor **Sediment** corresponds to three types of sediment (coded Sed1, Sed2, Sed3). For each type of Sediment, sampling was carried out on three sites (3 sites for Sed1, 3 sites for Sed2, 3 sites for Sed3).    * **Site** is coded : Site1, Site2, ..., Site9.   * The next factor is **Hydrology** : within each site, sampling is carried out in a dry plot and in a wet plot (coded Dry/Wet). Within each of the previous plot, sampling is carried out at two **Depths** (D1, D2) in triplicate. There are a total of n = 108 samples = 3 Sediment * 3 Sites * 2 Hydrology * 2 Depths * 3 replicates. I use the `lme()` function in R (nlme package) as follows :               Sediment <- as.factor(rep(c("Sed1","Sed2","Sed3"),each=36))     Site <- as.factor(rep(c("Site1","Site2","Site3","Site4","Site5",                                         "Site6","Site7","Site8","Site9"),each=12))     Hydrology <- as.factor(rep(rep(c("Dry","Wet"),each=6),9))     Depth <- as.factor(rep(rep(c("D1","D2"),each=3),18))     Variable <- rnorm(108)          mydata <- data.frame(Sediment,Site,Hydrology,Depth,Variable)          mod1 <- lme(Variable ~ Sediment*Hydrology*Depth, data=mydata,                   random=~1|Site/Hydrology/Depth)     anova(mod1)      I would like to run a _post-hoc_ comparison to test whether a term is significant or not. I'm able to do it for a simple main effect (e.g., **Sediment** ):               summary(glht(mod1,linfct=mcp(Sediment="Tukey")))      But the `glht()` function doesn't work for interaction terms. I found that the following could work for a 2-way anova :               mod1 <- lme(Variable~Sediment*Hydrology, data=mydata,                  random=~1|Site/Hydrology)     mydata$SH <- interaction(mydata$Sediment, mydata$Hydrology)     mod2 <- lme(Variable ~ -1 + SH, data=mydata, random=~1|Site/Hydrology)     summary(glht(mod2, linfct=mcp(SH="Tukey")))      Is it possible to use the same approach in the case of a 3-way anova? Any help on the way to do _post-hoc_ comparison on interaction terms in this case would be much appreciated.


title: Turning MiniBatchKMeans into Fuzzy MiniBatchKMeans

text: I'm using `Scikit-Learn`, which has an implementation of `MiniBatchKMeans`. I'm very inexperienced with ML, so I'm wondering how (if there is a way) I can adapt `MiniBatchKMeans` to be a fuzzy clustering algorithm.


title: How can I optimise computational efficiency when fitting a complex model to a large data set repeatedly?

text: I am having performance issues using the `MCMCglmm` package in R to run a mixed effects model. The code looks like this:               MC1<-MCMCglmm(bull~1,random=~school,data=dt,family="categorical"     , prior=list(R=list(V=1,fix=1), G=list(G1=list(V=1, nu=0)))     , slice=T, nitt=iter, ,burnin=burn, verbose=F)      There are around 20,000 observations in the data and they are clustered in around 200 schools. I have dropped all unused variables from the dataframe and removed all other objects from memory, prior to running. The problem I have is that it takes a very long time to run, unless I reduce the iterations to an unacceptably small number. With 50,000 iterations, it takes 5 hours and I have many different models to run. So I would like to know if there are ways to speed up the code execution, or other packages I could use. I am using `MCMCglmm` because I want confidence intervals for the random effects. On the other hand, I was hoping to get a new PC later this year but with a little luck I may be able to bring that forward, so I have been wondering how to best spend a limited amount of money on new hardware - more RAM, faster CPU etc. From watching the task manager I don't believe RAM is the issue (it never gets above 50% of physical used), but the CPU usage doesn't get much above 50% either, which strikes me as odd. My current setup is a intel core i5 2.66GHz, 4GB RAM, 7200rpm HDD. Is it reasonable to just get the fastest CPU as possible, at the expense of additional RAM ? I also wondered about the effect of level 3 CPU cache size on statistical computing problems like this ? **Update:** Having asked on meta SO I have been advised to rephrase the question and post on Superuser. In order to do so I need to give more details about what is going on "under the hood" in MCMCglmm. Am I right in thinking that the bulk of the computations time is spent doing optimisation - I mean finding the maximum of some complicated function ? Is matrix inversion and/or other linear algebra operations also a common operation that could be causing bottlenecks ? Any other information I could give to the Superuser community would be most gratefully received.


title: Are tree map diagrams effective at conveying information?

text: I have been looking at Treemapping as a potential data visualisation technique. Ref: Wikipedia article on Treemapping While popular in many visualisations, they don't seem to me to be effective at conveying information. This is because it seems to require attentive processing to understand and is difficult to make meaningful comparisons between elements. The mess seems to get in the way of a point. I'd like to know your thoughts: **Are treemaps effective at conveying information?**


title: Averaging binomial confidence intervals

text: I have a general statistical question. Three schools A, B, and C have 0.23, 0.56, and 0.8 fraction of female teachers. The upper and lower confidence intervals (CI) for the three schools can be tabulated as shown below:               school  female_instructors  upper_CI  lower_CI     A       0.33                0.53      0.13     B       0.55                0.65      0.45     C       0.80                0.90      0.70      If I were asked to report the average fraction of female instructors in the three schools combined I would have informed: (0.33+0.55+0.8)/3 or 0.56. How would I define the CI for the average fraction?? Average given CIs as well i.e. upper: 0.69 and lower:0.43?? Is this is the best way to do such calculations?


title: When is it appropriate to pool data?

text: Let's say I have some data on the amount of money spent on tv ads and total revenue from all sales. The data is available for three separate month.                                      Jan     Feb    March      $ spent on tv ads      100     110    150     total sales revenue    1000    1000   1500      The problem I am trying to solve is whether the marketing campaign (tv ads) had a significant effect on total sales revenue. Now, I'm still not sure on how to proceed in terms of tests, either time series analysis, etc. (suggestions?) But one other issue is going to occur, and that regards whether or not to pool the data. Should I run separate tests for each month or do I pool the data together? So my question is: When is it appropriate to pool data for analysis?


title: Hidden Markov model thresholding

text: I have developed a proof of concept system for sound recognition using mfcc and hidden markov models. It gives promising results when I test the system on known sounds. Although the system, when an unknown sound is inputted returns the result with the closest match and the score is not that distinct to devise it is an unknown sound e.g.: I have trained 3 hidden markov models one for speech, one for water coming out of water tap and one for knocking on the desk. Then I test them on unseen data and get following results:               input: speech     HMM\knocking:  -1213.8911146444477     HMM\speech:  -617.8735676792728     HMM\watertap:  -1504.4735097322673          So highest score speech which is correct          input: watertap     HMM\knocking:  -3715.7246152783955     HMM\speech:  -4302.67960438553     HMM\watertap:  -1965.6149147201534          So highest score watertap which is correct          input: knocking     HMM\filler  -806.7248912250212     HMM\knocking:  -756.4428782636676     HMM\speech:  -1201.686687761133     HMM\watertap:  -3025.181144273698          So highest score knocking which is correct          input: unknown     HMM\knocking:  -4369.1702184688975     HMM\speech:  -5090.37122832872     HMM\watertap:  -7717.501505674925     Here the input is an unknown sound but it still returns the closest match as there is no system for thresholding/garbage filtering.      I know that in keyword spotting an OOV (out of vocabulary) sound can be filtered out using a garbage or filler model but it says it is trained using a finite set of unknown words where this can't be applied to my system as I don't know all the sounds that the system may record. How is a similar problem solved in speech recognition system? And how can I solve my problem to avoid false positives?


title: How to perform a meta-analysis of studies without a control group?

text: I want to meta-analyse some studies which I've already collected and review. However, the studies did not used control groups, due to the types of questions used. How can I do the meta-analysis or which software can I use? Briefly my problem is as follows: I reviewed studies which the goal was to evaluate the perception of certain samples of their own behaviour. For each study, the proportion of the sample which claims to adopt the behaviour in analysis, is equivalent to the effect size. Do you have any idea how can I determine the average effect size? Thank you


title: Problem with Pareto distribution and R

text: I am trying to test this property of pareto distribution: Let f(x) be a pareto distribution $$ f(x)=\alpha \frac{x_m^\alpha}{x^{\alpha+1}} $$ so we have the cdf that is $$ CDF(x)=\int_{x_m}^{x}\alpha \frac{t_m^\alpha}{t^{\alpha+1}}dt=1-\frac{x_m^\alpha}{x^\alpha} $$ then the probability that $x>x_0$ is $$ P(x>x_0)=1-CDF(x)=\frac{x_m^\alpha}{x^\alpha} $$ and so we have $$ \frac{P(x>x_0)}{f(x)}=\frac{x}{\alpha} $$ Now i am trying to test it with R.                library(PtProcess)      dd<-rpareto(10000,1.5,0.01)      cdf<-ecdf(dd)      df<-density(dd)      ff<-(1-cdf(df$x))/df$y      If i plot ff                plot(df$x,ff)      I do not obtain the correct straight line. I guess that this is due at the way density() and ecdf() works. I need this form of the test (an a posteriori evaluation of fd and cdf) in order to perform the same test on a sample of data of unknown orgin. I guess that i need a way to binning the ecdf() function in the same way as hist() is the binning version of density. So my question is:   * Does there exist an equivalent binned function of ecdf() as hist() is the binned function of density()?    * or can I simulate ecdf() with hist()?


title: Proving no unbiased estimator exists for $\theta^{-1}$ for Poisson Dist. Why is $X_1,...,X_n = \sum{X_i}$?

text: Suppose that $X_1,...,X_n$ is an iid random sample from a Poisson distribution with mean $\theta$. I would like to prove that there exists no unbiased estimator of $\frac{1}{\theta}$. To do so, I will let $\delta(X)$ be an estimator of $\frac{1}{\theta}$. Then, I'd like to equate the expectation of $\delta(X)$ and $\frac{1}{\theta}$ E$[\delta(X)]$ = sum from x=0 to infinity of [ $\delta(x)$$P(X=x)$ ]. Now, my problem is that some books let $Y= \sum{X_i}$ and then have: E$[\delta(Y)]$ = sum from y=0 to infinity of [ $\delta(y)$$P(Y=y)$ ]. How are these two methods equivalent? How is it that the sum of the X's some how is equivalent to $X_1,...,X_n$? In other words, why is: $\delta(X_1,...,X_n) = \delta(Y)=\delta(\sum{X_i})$? Thanks!!!


title: R. function hclust, euclidean distances and cosines

text: My data is a table of cosines and I want to analyze it with hclust, which works on squared Euclidean distances. shall I do: d <- dist(mydata, method = "euclidean") fit <- hclust(d, method="ward") or rather: d=(dist(mydata))^2   fit <- hclust(d, method="ward") I don't know how to treat my cosines with hclust, which needs squared Euclidean distances. thanks!!!


title: slowing down effect

text: I am looking for an algorithm or function to look at the difference in reaction times after negative feedback. The problem is that I need to consider just the times in which negative feedback came after 3 positive feedbacks and no all the participants have the same number of negative feedback (or losses) Any idea? Thanks


title: Shapiro-Francia test error

text: I'm trying to run a normality test on the residuals after fitting a mixed- effect model (with lmer). I read that the Shapiro-Francia test can deal with data with more than 5000 observations (I have more than 8000), but when I run it I get an error:               sf.test(resid(dat.lmer11))      Error in sf.test(resid(dat.lmer11)) : sample size must be between 5 and 5000 Could anyone help, please?


title: How to implement rounding in BUGS in a survey context?

text: In a traditional filtering context, you'd analyse polls coming in with something like               mu[i] <- alpha[date[i]]     prec[i]<- SampleSize[i]/(alpha[date[i]]*(1-alpha[date[i]])     y[i] ~ dnorm(mu[i],prec[i])      If the surveys are observed with rounding, this _usually_ won't make a difference. But I'm in a situation where it does. What would be a good way to handle the observations being rounded in BUGS code?


title: Heteroskedasticity in linear regression model & data transformation

text: pls correct me if i'm wrong. of the econometrics' literature i've read so far, most mentioned heteroskedasticity is not a major problem empirically but multicollinearity would pose a greater concern to researchers and that data transformation would improve empirical results but not completely remove heteroskedasticity. Is this right? by above, i've used R to conduct a basic linear regression with one regressor and three dummies and conducted Breusch-Pagan test using package 'car' and VIF test using package 'lmtest'. Thereafter, i conducted a second regression, log- transforming the dependent and independent variable and carried similar tests as aforementioned. for both models, m1 for former and m2 for latter, and using VIF test, both presented less than 5 and reasonably to say multicollinearity isn't a problem however, for Breusch-Pagan, m1 failed to reject null but m2 (logged) rejected it. can someone enlighten me why the result is not consistent with the theory above? some related questions,   1. would it better not to log-transform the model?   2. i'm suspecting that the regressor is correlated with the residuals. is there test(s) 'out there' for endogeneity?   3. if log model is preferred, what can i then do to 'lessen' heteroskedasticity in m2? Thanks in advance! Dylan


title: Start time requirements or assumptions for survival analysis

text: We have prospective data from an observational registry and wish to consider the affects of a gene on time to cardiovascular events. The data includes standard data like age, gender, ... and also the date(age) of a cardio event. In general terms what are the conditions/assumptions/requirements of a starting time for survival analysis? Entry into the registry is at random times for each individual (as he or she is recruited into the registry). This is obviously no data at which the registry participants' genes changed, and there no intervention. Is it valid to use the date of enrollment in the registry as the start time for the survival analysis? We (group working on this) are asking as some work has been done and the the enrollment date was used as the start time for the survival interval and there is some disagreement now as to whether this is a valid approach.


title: Does corr function in Matlab consider correction for tied ranks?

text: Does 'corr' function in Matlab consider correction for tied ranks while calculating Spearman correlation coefficient? Thanks.


title: Probability of picking a biased coin

text: Suppose you have a bag of 100 coins of which 1 is biased with both sides as Heads. You pick a coin from the bag and toss it three times. The result of all three tosses is Heads. What is the probability that the selected coin is biased? My answer:-               P(selecting a biased coin) = 1/100     P(getting a head thrice with the biased coin) = 1          P(selecting an unbiased coin) = 99/100     P(getting a head thrice with the unbiased coin) = 1/8          P(selecting a biased coin|coin toss resulted in 3 heads) =                P(selecting a biased coin and getting heads thrice)/               [P(selecting a biased coin and getting heads thrice) +                P(selecting an unbiased coin and getting heads thrice)]           = (1/100)/[(1/100) + (99/800)]     = (1/100)/(107/800)     = 8/107     = 0.0747      Is this correct? Thanks.


title: About canonical correlation analysis

text: > Given two column vectors $X = (x_1, \dots, x_n)'$ and $Y = (y_1, \dots, > y_m)'$ of random variables with finite second moments, canonical-correlation > analysis seeks vectors $a$ and $b$ such that the random variables $a' X$ and > $b' Y$ maximize the correlation $\rho = \operatorname{corr}(a' X, b' Y) = > \displaystyle\frac{a' \Sigma _{XY} b}{\sqrt{a' \Sigma _{XX} a} \sqrt{b' > \Sigma _{YY} b}}$. > > (Wikipedia) Now, at times, I see this written into the form > $max \ a' \Sigma _{XY} b \\ s.t.\ a' \Sigma_{XX} a = b' \Sigma_{YY} b = 1$ I do see that scaling $a$ or $b$ with any scalar has no effect in the first equation, but how do I get to the second equation?


title: Get distribution curve characteristics by selecting peak point?

text: I have no math or stats background, so this may be a real dumbo question.   I'm developing a time estimation app which would let user enter her worst case and best case estimates, and then select a probability curve which she feels looks realistic. I needed to draw some kind of curve that sums up to 1, so I came up with a normal distribution. I'm not sure that it captures how we give time estimates but this is not what question is about. I want user to be able to move finger on the chart, so the curve would bend accordingly. I assume it would feel nice if curve's peak followed the finger, but then I have the following question: is it possible to get normal distribution curve by arbitrary peak point? If not, which curve is good for this? **UPDATE:**   I see that I can obtain _μ_ from _x_ location.   Can I get _σ_ for an arbitrary _y_ though (combined with known _μ_ )?


title: Correlation of two 1d point processes with non-uniform unconditional distribution

text: I want to measure the correlation between two 1D point processes $x$ and $y$. Ordinarily I could use the bivariate K-function $K(t) = \frac{T}{n_xn_y} \sum_{i=1}^{n_x} \sum_{j=1}^{n_y} w(x_i,y_j) I[d(x_i,y_j)<t]$ where $n_x$ is the number of observations in $x$ and $n_y$ is the number of observations in $y$. Deviation from $K(t)=t$ is an indication of correlation between the two point processes. However, my analysis is complicated by the fact that the unconditional distributions of the point processes are non-uniform. They each have a characteristic distribution of intervals between observations, which could be different for $x$ and $y$. To make things more complicated, the intensities vary over the time period (for example, the intensity might be low around midnight and higher during daylight hours). Assuming I have enough data to get good estimates of the unconditional distributions of $x$ and $y$, and of the way that intensities vary across time, is there a way to take this into account?


title: Predicting long-memory processes

text: I'm working with a two-state process with $x_t$ in $\{1, -1\}$ for $t = 1, 2, \ldots$ The autocorrelation function is indicative of a process with long-memory, i.e. it displays a power law decay with an exponent < 1. You can simulate a similar series in R with:               > library(fArma)     > x<-fgnSim(10000,H=0.8)     > x<-sign(x)     > acf(x)      My question: is there a canonical way to optimally predict the next value in the series given just the autocorrelation function? One way to predict is simply to use $\hat{x}(t) = x(t-1)$ which has a classification rate of $(1 + \rho_1) / 2$, where $\rho$ is the lag-1 autocorrelation, but I feel like it must be possible to do better by taking into account the long-memory structure.


title: Alternative ways for interpretation of odds

text: I am reading the second edition of Categorical Data Analysis by Alan Agresti, and somehow stuck in the following second paragraph: ![enter image description here](http://i.stack.imgur.com/ahQDT.png) I don't quite understand why $\beta\pi(\hat{x})(1 - \pi(\hat{x}))$ will give the probability when $x = 26.3$, can anyone enlighten me? Thanks.


title: Generic test for unimodality given sample

text: Are there any generic tests to validate if a given sample follows a unimodal distribution, like a Gaussian, Cauchy, Student's t or a chi-square?


title: How to set confidence level for wilcoxsign_test (package coin)?

text: In R, the function `wilcox.test` takes the argument `conf.level = 0.095` (for example). Giving the same argument to the function `wilcoxsign_test` from the _coin_ package returns a warning:               additional arguments conf.level will be ignored      **What default confidence level does`wilcoxsign_test` use, and how can I change it?** _Or :_ Why do I not need a confidence level for this function?


title: Limit of quantile normalization

text: Is quantile normalization adequate for normalizing data with very few samples? For example this microarray data. Typically after normalization we'd like to compare `Cancer-1 with Normal` and `Cancer-2 with Normal` for differential expression.                   mRNA     Cancer-Type-1   Cancer-Type-2  Normal         -----------------------------------------------------------         mRNA1      30        49    12         mRNA2     199        200   78         ...        ...       ...  ....         mRNA1000   13        40    88      If not what is the appropriate normalization method?


title: How would re-weighting American Community Survey diversity data affect its margins of error?

text: Background: My organization currently compares its workforce diversity statistics (ex. % persons with disabilities, % women, % veterans) to the total labor force availability for those groups based on the American Community Survey (a surveying project by the US Census Bureau). This is an inaccurate benchmark, because we have a very specific set of jobs that have different demographics than the labor force as a whole. Say, for example, that my organization is mostly engineers. Engineering is only about 20% women in my state. If we compare ourselves to the total labor force benchmark, which is more like of 50% women, it results in panic that “we only have 20% women, this is a disaster!” when really, 20% is what we should be expecting because that’s what the labor landscape looks like. My goal: What I would like to do is take the American Community Survey occupation data (by diversity category) and re-weight it based on the composition of jobs in my business. Here is a sample data set for Social and Community Service workers. I want to add these job codes listed together (because our crosswalk is to job groups, not to specific job codes), then I want to weight that benchmark based on the number of people we have in that category (ex. our 3,000 Social and Community Service workers), then I want to do the same to all the other job groups, add those numbers together, and divide by our total number of workers. This would give me a new re-weighted diversity measure (ex. from 6% persons with a disability to 2% persons with a disability). My questions: How do I fit margins of error to this final rolled up benchmark? I do not have the raw census data set (obviously), but you can view margins of error for each number in the link that I provided by toggling the “Estimate” field to "Margin of Error" at the top of the table. My other co-workers who are working with this data fully intend to ignore the margins of error, but I am worried that we are creating a statistically meaningless benchmark for ourselves. Is this data even still usable after the manipulation described above?


title: How to report significance of a variable that is dependent on two others?

text: Cost per click is defined by the cost per impression and the click through rate. Given a period of time, divide number of clicks by the sum of the cost per impression for all ads in the same time period. Cost per click is equivalent to click through rate multiplied by cost per impression (which are assumed to be independent). I already know how to report significance for changes in click through rate and for cost per impression. How can I report a significance value for cost per click?               cpi = cost / impressions      ctr = clicks / impressions     cpc = cost / clicks          cpc / clicks = cost     cpc / clicks / impressions = cost / impressions     cpc = (clicks / impressions) * (cost / impressions)     cpc = ctr * cpi


title: Variable Selection One by One vs Simultaneously

text: The high dimensional variable selection problem is really popular now. But I have a question: If I do simple linear regression regressing one response variable on 1 covariate at a time first and then control the FDR to select the significant feature variables, what's the disadvantage of this comparing to the lasso or group lasso algorithm, which choose the feature variables simultaneously? Basically the question could be reduced to **what's the difference between regressing one response variable on multiple covariates vs on 1 covariate at a time?**


title: How do I perform a statistical analysis on an intervention based study?

text: I am trying to get my head around the best way to analyse a study. The study comprises of persons that get given a number of up to 10 different interventions, whereby hospital admission frequency is measured as an outcome. Each intervention are not given out randomly, and some persons are not required to have all 10. Also, the interventions can be given out in different orders for each person, but a given intervention may only be used once for each person. **As such, it is a given that a majority of these people will be living with, and experiencing the effects of multiple interventions at the same time** The questions is then to measure the effect of each individual intervention by looking at the hospital outcome data? I am finding this a tricky concept to get my head around as for example if I want to measure the effect of intervention 1, I might say what is the total amount of person time where people are living before intervention, and the total amount of person time people are living after the intervention - where I measure the health outcome of both periods. The problem is that some of these persons years will also have had persons years living before/after a different intervention. Does anybody have any ideas on how I can use the data to measure the effect of each individual variable?


title: Parametric vs. Nonparametric

text: I did perform a search for this question (I thought it was bound to be asked), but I haven't found one; hopefully this won't be a duplicate. I'm trying to decide if I should take a course in stats that is more parametric or one that is nonparametric. From a maths perspective, what would be the advantages of learning one over the other? Would one be more advantageous over the other as far as working in industry?


title: Clustering and correlations

text: I have a dataset of $n$ $p$-dimensional vectors (objects) that I want to cluster.   1. One way to do this is to compute the ($n \times n$) correlation matrix $C$, then obtain a dissimilarity matrix, $D$, from $C$ such that each element $d_{ij}$ in $D$ is a function of the single element $c_{ij}$ in $C$ (e.g., $D=1-C$), and then cluster on $D$.   2. Instead, I want to take $C$ and obtain $D$, such that each element $d_{ij}$ of $D$ is a function of the entire vectors $C_i$ and $C_j$ in $C$ (for example, $d_{ij}$ can be the Euclidean distance between $C_i$ and $C_j$). Why do 2)?:   * Approach 1) computes the dissimilarity (distance) between two input vectors $i$ and $j$ solely as a function of the similarity between $i$ and $j$; in contrast, approach 2) computes the dissimilarity (distance) between two input vectors $i$ and $j$ as a function of the similarity of the similarities between $i$ and all other vectors vs. $j$ and all other vectors.   * On the problems with which I am working, approach 2) seems to perform better in practice. What I am wondering **is whether there is any reason why** approach 2) would not be valid?


title: How to correlate graph properties with vertex attributes?

text: I have a directed graph, **S** , where each vertex stores several integers of attributes. I think it might be called a "attribute relational graph" or just an old fashioned "attributed graph". The networkx Python package allows any attributes to be added to each vertex, but I can't figure out how to pair these attributes with features of the graph. Let's say I have the following vertices representing moments in time:   * **S1** , where Alice has \$1, Bob has \$2, and Cindy has \$5.   * **S2** , where Alice has \$0, Bob has \$3, and Cindy has \$5.   * **S3** , where Alice has \$0, Bob has \$0, and Cindy has \$8. Edges between the vertices imply a temporal connection between the moments in time, for example, before and after a transaction between Alice, Bob, and Cindy. Additionally, vertex weight and edge weight are the frequency of that state or transaction occurring over a long period of sampling at discrete time intervals. My question is, how can I compute correlations between those numeric attributes and the graph topology? The only place I can think to start would be to cluster **S** , using the graph topology, and analyze the change of attributes between vertices within each cluster. For that you could define a similarity metric like the difference between the edge connecting each pair of vertices and check if this is correlated within a cluster. Are there any standard algorithms for this type of analysis or am I looking at all custom code? Other things I'd also be interested in computing,   * A greedy "loss-only" path for Alice, what is the shortest number of transactions to make her broke from any given point.   * The number of occurrences in the graph where Bob acts an an intermediate party to Alice and Cindy.


title: Probability of Two Samples Containing Overlapping Data

text: Quick question about sampling a data population that contains non-overlapping data: if I take two samples of 50 from a population of 100,000 unique data points, how would I calculate the probability of the two samples containing a shared data point? As well, how do I determine the population size required for the probability of overlap to be below a certain point (0.05%) for $n=50$? Thanks.


title: Simulate the distribution of the median

text: ![Question](http://i.stack.imgur.com/ujlEW.png) I have done part(a); now I need to do part(b): I want to simulate the distribution of the median $M$ using R. I don't know how to work it out. And for part(c), do I need to use R as well? This is what I have done               medians <- vector(length=100)     for (i in 1:100) {          x <- rnorm(n=15,mean=0,sd=1)          medians[i] <- median(x)      }          hist(medians)


title: Topic modeling for grouped data

text: Is there a variation of Latent Dirichlet Allocation (LDA) for grouped data? As an example, let us consider corpus of all Yahoo Q&A (where for simplicity we consider a question lumped together with all the corresponding answers as a document). While submitting a question to Yahoo's site, user has to specify a category to which it belongs to (let us, again for simplicity assume that categories are simple, instead of hierarchical). Is there a version of LDA that takes advantage of this information?


title: How to derive the conditional posterior density in hierarchical bayesian models?

text: I was reading on Gelman's Bayesian Data Analysis - Chapter 5 - Hierarchical model Suppose: data : $y_j$ s parameter: $\theta$ hyperparameter: $\phi$ On page 126, he mentions the analytical derivation of conditional and marginal distributions. I understand that we can easily derive the joint posterior distribution using factorisation based on Bayes' theorem, but I don't understand how to derive the conditional posterior distribution $p(\theta|\phi ,y)$. Is there a general method to go about deriving the conditional posterior distribution? Based on Bayes' theorem perhaps? Please explain thoroughly.


title: How to set up a posterior predictive test quantities (Bayesian context) to check for independent Poisson distributions?

text: Suppose we are given data $y_j \sim \text{Poi}(\lambda)$ and assume $y_j$ are iid. We can assume the prior distribution for $\theta$ follows $\text{Gamma}(\alpha, \beta)$, The posterior distribution of $\theta | y$ is also Gamma while the posterior predictive distribution is Negative Binomial. My question is that after deriving everything above, how would you go about setting up a posterior predictive **test quantities** to check for **_independent_** Poisson distributions? Notice that test quantity is a function of parameters and data in a Bayesian context. It is similar to a test statistic in classical statistic.


title: Machine learning techniques for parsing strings?

text: I have a lot of address strings:               1600 Pennsylvania Ave, Washington, DC 20500 USA      I want to parse them into their components:               street: 1600 Pennsylvania Ave     city: Washington     province: DC     postcode: 20500     country: USA      But of course the data is dirty: it comes from many countries in many languages, written in different ways, contains misspellings, is missing pieces, has extra junk, etc. Right now our approach is to use rules combined with fuzzy gazetteer matching, but we'd like to explore machine learning techniques. We have labeled training data for supervised learning. The question is, what sort of machine learning problem is this? It doesn't really seem to be clustering, or classification, or regression.... The closest I can come up with would be classifying each token, but then you really want to classify them all simultaneously, satisfying constraints like "there should be at most one country;" and really there are many ways to tokenize a string, and you want to try each one and pick the best.... I know there exists a thing called statistical parsing, but don't know anything about it. So: what machine learning techniques could I explore for parsing addresses?


title: Deriving Negentropy. Getting stuck

text: So, this question is somewhat involved but I have painstakingly tried to make it as straight-forward as possible. **Goal:** Long story short, there is a derivation of negentropy that does _not_ involve higher order cummulants, and I am trying to understand how it was derived. **Background: (I understand all this)** I am self-studying the book 'Independent Component Analysis', found here. (This question is from section 5.6, incase you have the book - 'Approximation of Entropy by Nonpolynomial Functions'). We have $x$, which is a random variable, and whose Negentropy we wish to estimate, from some observations we have. The PDF of $x$ is given by $p_x(\zeta)$. Negentropy is simply the difference between the differential entropy of a standardized gaussian random variable, and the differential entropy of $x$. The differential entropy here is given by $H$, such that: $$ H(x) = -\int_{-\infty}^{\infty} p_x(\zeta) \: log(p_x(\zeta)) \: d\zeta $$ and so, the negentropy is given by $$J(x) = H(v) - H(x)$$ where $v$ is a standardized gaussian r.v, with PDF given by $\phi(\zeta)$. Now, as part of this new method, my book has derived an estimate of the PDF of $x$, given by: $$ p_x(\zeta) = \phi(\zeta) [1 + \sum_{i} c_i \; F^{i}(\zeta)] $$ (Where $c_i = \mathbb{E}\{F^i(x)\}$. By the way, $i$ is _not_ a power, but an index instead). For now, I 'accept' this new PDF formula, and will ask about it another day. This is not my main issue. What he does now though, is plug this version of the PDF of $x$ back into the negentropy equation, and ends up with: $$ J(x) \approx \frac{1}{2}\sum_i\mathbb{E} \{F^i(x)\}^2 $$ Bear in mind, the sigma (here and for the rest of the post), just loops around the index $i$. For example, if we only had two functions, the signal would loop for $i=2$ and $i=2$. Of course, I should tell you about those functions he is using. So apparently, those functions $F^i$ are defined as thus: > The functions $F^i$ are not polynomial functions in this case. (We assume > that the r.v. $x$ is zero mean, and of unit variance). Now, let us make some > constraints and give properties of those functions: > > $$ F^{n+1}(\zeta) = \zeta, \: \: c_{n+1} = 0 $$ > > $$ F^{n+2}(\zeta) = \zeta^2, \: \: c_{n+1} = 1 $$ > > To simplify calculations, let us make another, purely technical assumption: > The functions $F^i, i = 1, ... n$, form an orthonormal system, as such: > > $$ \int \phi(\zeta) F^i(\zeta)F^j(\zeta)d\zeta= \begin{cases} 1, \quad > \text{if } i = j \\ 0, \quad \text{if } i \neq j \end{cases} $$ > > and > > $$ \int \phi(\zeta)F^i(\zeta)\zeta^k d(\zeta) = 0, \quad \text{for } k = > 0,1,2 $$ Almost there! Ok, so all that was the background, and now for the question. The task is to then, simply place this new PDF into the differential entropy formula, $H(x)$. If I understand this, I will understand the rest. Now, the book gives the derivation, (and I agree with it), but I get stuck towards the end, because I do not know/see how he is cancelling out. Also, I do not know how to interpret the small-o notation from the taylor expansion. This is the result: Using the taylor expansion $(1+\epsilon)log(1+\epsilon) = \epsilon + \frac{\epsilon^2}{2} + o(\epsilon^2)$, for $H(x)$ we get: $$ H(x) = -\int \phi(\zeta) \; (1 + \sum c_i F^i(\zeta)) \; (log(1 + \sum c_i F^i(\zeta) + log(\zeta)) \; d(\zeta) \\ = -\int \phi(\zeta) log(\zeta) -\int \phi(\zeta) \sum c_i F^i(\zeta) log(\phi(\zeta)) -\int \phi(\zeta) \; [\sum c_i F^i(\zeta) + \frac{1}{2}(\sum c_i F^i(\zeta))^2 + o((\sum c_i F^i(\zeta))^2)] $$ and so **The Question: (I dont understand this)** $$ H(x) = H(v) - 0 - 0 -\frac{1}{2}\sum c_i^2 + o((\sum c_i)^2 $$ So, my problem: Except for the $H(v)$, I dont understand how he got the final 4 terms in the last equation. (ie, the 0, the 0, and the last 2 terms). I understand everything before that. He says he has exploited the orthogonality relationships given in the properties above, but I dont see how. (I also dont understand the small-o notation here, in the sense of, how it is used?) THANKS!!!! **EDIT:** I have gone ahead and added the images from the book I am reading, it pretty much says what I said above, but just incase someone needs additional context. ![enter image description here](http://i.stack.imgur.com/WyAPQ.jpg) ![enter image description here](http://i.stack.imgur.com/4ZaFH.jpg) ![enter image description here](http://i.stack.imgur.com/STPx9.jpg) And here, marked in red, is the exact part that is confusing me. How does he use the orthogonality properties to get that last part, where things are cancelling out, AND the final summations involving $c_i^2$, and the small-o notation summation? Thanks.


title: How to apply published factor analysis loadings to sample data to compute factor scores?

text: I have data on 18 participants responding to a 30-item instrument. I also have a paper in which the researchers performed a principal axis factor analysis on the same survey where they report the item loadings (4 factors). **Is it possible to convert the data that I have, and compute scores for each participant on the four factors based on the published loadings?**


title: Determining if awards are fair

text: Imagine a population of people where $x\%$ have blue hats and $(100-x)\%$ have green hats where every year some people are awarded prizes. The rule is that everyone should have an equal chance of being awarded a prize every year. The village bosses would like to decide if this rule is being followed as they suspect there may be some prejudice based on hat color. They count the number of people with blue hats who get a prize and call it $A$ and also the number of people with green hats who get a prize and call that $B$. The total population size of green and blue hats is $P$. What is a good method to determine if the prize awarding has indeed been fair?


title: Non-Correlated errors from Generalized Least Square model (GLS)

text: As a financial institution, we often run into analysis of time series data. A lot of times we end up doing regression using time series variables. As this happens, we often encounter residuals with time series structure that violates basic assumption of independent errors in OLS regression. Recently we are building another model in which I believe we have regression with autocorrelated errors.The residuals from linear model have `lm(object)` has clearly a AR(1) structure, as evident from ACF and PACF. I took two different approaches, the first one was obviously to fit the model using Generalized least squares `gls()` in R. My expectation was that the residuals from gls(object) would be a white noise (independent errors). But the residuals from `gls(object)` still have the same ARIMA structure as in the ordinary regression. Unfortunately there is something wrong in what I am doing that I could not figure out. Hence I decided to manually adjust the regression coefficients from the linear model (OLS estimates). Surprisingly that seems to be working when I plotted the residuals from adjusted regression (the residuals are white noise). I really want to use `gls()` in `nlme` package so that coding will be lot simpler and easier. What would be the approach I should take here? Am I supposed to use REML? or is my expectation of non- correlated residuals (white noise) from gls() object wrong?               gls.bk_ai <- gls(PRNP_BK_actINV ~ PRM_BK_INV_ENDING + NPRM_BK_INV_ENDING,                       correlation=corARMA(p=1), method='ML',  data  = fit.cap01A)          gls2.bk_ai  <- update(gls.bk_ai, correlation = corARMA(p=2))          gls3.bk_ai <- update(gls.bk_ai, correlation = corARMA(p=3))          gls0.bk_ai <- update(gls.bk_ai, correlation = NULL)          anova(gls.bk_ai, gls2.bk_ai, gls3.bk_ai, gls0.bk_ai)            ## looking at the AIC value, gls model with AR(1) will be the best bet          acf2(residuals(gls.bk_ai)) # residuals are not white noise      Is there something wrong with what I am doing???????


title: Confusion related to density estimation

text: I have this confusion related to how this density was estimate. I have attached the screenshot of the paper Any references related to basis expansion and all. I didn't get how this was derived actually. ![enter image description here](http://i.stack.imgur.com/R5Xxu.png)


title: How to include an interaction term in GAM?

text: The following code evaluates the similarity between two time series:               set.seed(10)     RandData <- rnorm(8760*2)     America <- rep(c('NewYork','Miami'),each=8760)          Date = seq(from=as.POSIXct("1991-01-01 00:00"),                 to=as.POSIXct("1991-12-31 23:00"), length=8760)          DatNew <- data.frame(Loc = America,                          Doy = as.numeric(format(Date,format = "%j")),                          Tod = as.numeric(format(Date,format = "%H")),                          Temp = RandData,                          DecTime = rep(seq(1, length(RandData)/2) / (length(RandData)/2),                                        2))     require(mgcv)     mod1 <- gam(Temp ~ Loc + s(Doy) + s(Doy,by = Loc) +       s(Tod) + s(Tod,by = Loc),data = DatNew, method = "ML")      Here, `gam` is used to evaluate how the temperature at New York and Miami vary from the mean temperature (of both locations) at different times of the day. The problem that I now have is that I need to include an interaction term which shows how the temperature of each location varies throughout at the day for different days of the year. I eventually hope to display all of this information on one graph (for each location). So, for Miami I hope to have one graph that shows how the temperature varies from the mean during different times of the day and different times of the year (3d plot?)


title: Confusion related to expectation

text: I was reading a paper related to minimaxity. ![enter image description here](http://i.stack.imgur.com/5tTGZ.png) If $\hat{\theta}$ is the estimate of an unknown parameter $\theta$ and $l(\hat\theta-\theta)$ is the loss function they have mentioned $E_p$ is the mean of $l(\hat\theta-\theta)$ under distribution P. I didn't get this part. If $\hat\theta$ is an estimate and $\theta$ is the unknown parameter of a distribution P, then $\theta$ is fixed since it is the mean of the distribution P which is not known. So what is meant by expectation of $l(\hat\theta-\theta)$


title: How would I find the P-Value of this data?

text: I have the following data (relating to consecutive months at current job):   * Mean: 45.4   * Standard Deviation: 60.89   * Sample Size: 48 I am asked to do a one-tailed test at a 1% significance level and provide:   1. the value of the test statistic (from above).   2. P-Value of the test statistic. My initial assumption about the data is that most of the values would fall under 60 months. With that said, my hypothesis is: H0: x̄ ≤ 60, Ha: x̄ > 60 How do I proceed with this problem? I am not sure I understand the process. Thank you


title: Maximal & closed frequent -- Answer Included

text: $$My \ \ dataset:$$ $$1: A,B,C,E$$ $$2:A,C,D,E$$ $$3:\ \ \ \ \ B,C,E$$ $$4:A,C,D,E$$ $$5:\ \ \ \ C, D, E$$ $$6: \ \ \ \ A, D,E$$ I want to find out the **maximal frequent item sets** and the **closed frequent item sets**.   * Frequent item set $X ∈ F$ is **maximal** if it does not have any frequent supersets.   * Frequent item set X ∈ F is **closed** if it has no superset with the same frequency * * * So I counted the occurrence of each item set.               {A} = 4 ;  {B} = 2  ; {C} = 5  ; {D} = 4  ; {E} = 6          {A,B} = 1; {A,C} = 3; {A,D} = 3; {A,E} = 4; {B,C} = 2;      {B,D} = 0; {B,E} = 2; {C,D} = 3; {C,E} = 5; {D,E} = 3          {A,B,C} = 1; {A,B,D} = 0; {A,B,E} = 1; {A,C,D} = 2; {A,C,E} = 3;      {A,D,E} = 3; {B,C,D} = 0; {B,C,E} = 2; {C,D,E} = 3          {A,B,C,D} = 0; {A,B,C,E} = 1; {B,C,D,E} = 0      * * * Min_Support set to $50%$ // Very important. Thanks steffen for reminding of that. Does **maximal** = $\{{A,B,C,E\}}$ ? Does **closed** = $\{{A,B,C,D\}} \ and \ \{{B,C,D,E\}}$?


title: Independent or dependent data

text: I have lots of football result data, and I want to test the dependency of the number of goals per game on the weather (e.g. sun, rain, snow). I'm struggling to decide if the data is dependent or independent (and hence which ANOVA test to use). I can't be sure that every team will have played every other team more than once.


title: Decision rules based on amount of overlap between Bayesian probability intervals

text: My research hypothesis is that an intervention, I, results in greater scores than the control,C. If the intervals for these two don't overlap, I know I can make a clear-cut statement either accepting or rejecting my hypothesis. However, I also want to make "degree of evidence in favor of" types of statements should they overlap. Below would be my decision rule: Assuming the width of the intervals is similar enough to not bias the comparison, if the upper limit for I is greater than the upper limit for C and the overlap, O, is 0% < O < 50%, then let the research hypothesis be supported yet penalized by O. If the above is true except 50% <= O < 100%, then there is too much evidence against my hypothesis to accept even with penalty. In all other cases of overlap my hypothesis would be rejected. This is intuitive to me but wasn't to my stats professor. It's not that she disagreed insomuch as she was unsure. I like it because it seems similar to testing against a null model with an alternative with Bayes Factors where you have this relative support for and against, except I'm estimating parameters directly. What are some thoughts on this idea?


title: Multiple ANOVAs, shall I correct the homogeneity of variance test?

text: I want to run nine different ANOVAs on different variables that may change among groups. The ANOVAs are preceded by the Levene test of variance homogeneity. Of course, I want to correct for multiple comparisons in the final alpha threshold for the ANOVAs. But as random p values can be under 0.05 in all tests, the same should happen to the Levene test, am I right? So my question is, shall I correct the Levene test results also for multiple comparisons? What argument can it be not to correct the homogeneity test the same was the results from those data are corrected for multiple comparisons?


title: Bayesian analysis of an opinion poll

text: I'm sorry in advance for rather long questions. This is an example in "Bayesian logical data analysis for physical sciences" by P. C. Gregory and I have some questions about the example. > In a poll of 800 decided voters, 440 voters supported the political party A. > Let's denote the poll result as $D$. The quantity of interest is the > probability that the party A will achieve a majority of at least 51% in the > upcoming election, assuming the poll will be representative of the > population at the time of the election. The book regards the problem as a model selection problem. $M_1$ : The party A will achieve a majority with a parameter $H$ that has uniform prior in the range $0.51 \le H \le 1$. $M_2$ : The party A will not achieve a majority with a parameter $H$ that has uniform prior in the range $0 \le H < 0.51$. If we have no prior reason to prefer $M_1$ over $M_2$, we can write the odds ratio: $$ \begin{aligned} O_{12}&=p(M_1|D,I)/p(M_2|D,I)\\ &=p(D|M_1,I)/p(D|M_2,I)\\ &=\frac{\int_{0.51}^1 p(H|M_1,I)p(D|H,M_1,I) dH }{\int_{0}^{0.51} p(H|M_2,I)p(D|H,M_2,I) dH}\\ &=\frac{\int_{0.51}^1 (1/0.49)p(D|H,M_1,I) dH }{\int_{0}^{0.51} (1/0.51)p(D|H,M_2,I) dH}\\ &=87.68 \end{aligned} $$ Here are my questions: The book doesn't give explicit expressions for $p(D|H,M_1,I)$ and $p(D|H,M_2,I)$. If I use binomial distribution: $$ p(D|H,M_1,I)=p(D|H,M_2,I)=\frac{800! H^{440}(1-H)^{800-440}}{440!(800-440)!} $$ I get $87.03$ as a result. It is not same to the value $87.68$ of the Book. **What probability distribution should I use for the likelihoods?** **Why do I have to introduce the models $M_1$ and $M_2$? Why is is this not appropriate to approach the problem like this:** $$ O_{12}=\frac{\int_{0.51}^1 p(H|D,I) dH}{\int_{0}^{0.51} p(H|D,I) dH} $$ It does not have the factor $(1/0.49)/(1/0.51)$ introduced with the models $M_1$ and $M_2$.


title: Two-sided p-value of Fisher's exact test

text: For a paper I'm using Fisher's Exact Test in R. As I already found out the p.value in R is calculated as the sum of all probabilities of the hypergeometric distribution with the given parameters lower or equal the the probability of the observed value. And I already learned that there are multiple ways of calculating the two-sided p-value and that there is a ongoing discussion about this. So now my question is, why is R doing it this way? Is there any source I can refer to in my paper? FYI: In R the computation looks like this:               set.seed(42)     A <- c(sample(c(1,2),replace=T,size=15),1,1,1,1,1,1)     B <- c(sample(c(1,2),replace=T,size=15),1,1,1,1,1,1)     (tab <- table(A,B))     fisher.test(tab)$p.value          n11 <- tab[1,1]     white <- sum(tab[1,])     black <- sum(tab[2,])     drawn <- sum(tab[,1])          lim <- dhyper(n11,m=white,n=black,k=drawn)     all <- dhyper(1:drawn,m=white,n=black,k=drawn)     sum(all[all<=lim])


title: Baseline differences in RCT: Which variables (if any) should be included as covariates?

text: I recently completed a study whereby I randomly assigned participants to one of two treatment groups. I tested participants at baseline, immediately post- intervention, 1 months, and 4 months on a somewhat large number of outcome variables. I was planning on running several mixed ANOVAs to examine group x time interactions. Some of the comparisons will be 2 (group) x 2 (time: baseline and post-intervention) comparisons and some will be 2 (group) x 3 (time: baseline, 1 month, 4 month) comparisons. Before beginning my analyses, I compared the two treatment groups on all baseline variables. I found that the groups differ on 4 baseline variables if I use an alpha level of .05 or 2 baseline variables if I use an alpha level of .01 to compare the groups. I have two questions about this:   1. What alpha level should I be using to compare the groups at baseline? I was thinking an alpha level of .01 because I am comparing the two groups on 24 baseline characteristics and I thought I should chose a more stringent alpha level than .05 to reduce family-wise error rate seeing as a large number of tests are being performed, but from my readings it seems most people use .05. What do you recommend?   2. What do I do about these differences? I could include these variables as covariates, but my sample size is quite small and using 4 covariates does not seem appropriate (which is also partly why I am favouring only accepting differences if they are significant at the .05 level) Any help on this would be very much appreciated!


title: what to do with 0.5 class probabilities ?

text: I am currently training a random forest regressor (scikit learn) on the Titanic dataset. My question is related to this issue (http://stackoverflow.com/questions/19984957/scikit-predict-default-threshold) on stack overflow. I noted that I didn't have the same value as in scikit for measures like Precision, Recall, F1-score ... After investigating I noticed that the reason was I considered 0.5 probabilities individuals to be in class 1 while scikit classes them as 0. So here are my questions :   * is it better to class 0.5 probabilities individuals in 0 or 1 class ? On titanic for example it can change significantly the value of such measures.   * would it be legit not to use these ? I do not think so because it bias your results and may tend to improve them.    * what about classification with more than two classes ? If I have 1/3,1/3,1/3 as probabilities for one individual what should I do ?    * is there any performance measure emancipated from this problem ?    * is scikit-learn choosing this 0.5 -> 0 class every time or can it be random / depends on the model selected ?


title: Error when using corAR1 and varExp in mutilevel model

text: I'm using nlme in R to fit a multilevel model for some physiological (skin conductance level) responses while watching a film. I'm specifying the model as: model <- lme(SCL <- variableA * variableB + TIME, random=~TIME|SUBJECTS, data=data5, na.action=na.omit, control=list(opt="optim"), method="ML", correlation=corAR1(0.9)) So far, so good. But, since the error structure shows heterocedasticity, I also included the argument: weights=varExp(form=~TIME|SUBJECTS) The problem is that when I include this function (weights), I get the following error message: Error in recalc.corAR1(object[[i]], conLin): Na/NaN/Inf en llamada a una funcion externa (arg 1) I've looked in the help forums, but I haven't been able to find a solution. According to Bliese (Multilevel modeling in R), you can use together the corAR1 function (for correcting autocorrelation) and varExp (for heteroscedasticity), and it is Ok when I use them separately, but when I use both in a model I obatin the error above. Could anyone help me with this?? Thank you so much!


title: Estimating the size of an intersection of multiple sets by using a sample of one set

text: I'm working on an algorithm which needs to calculate the size of a set generated by the intersections of at least 2 sets. More specifically: $$ z = \left |A_0 \cap \ldots \cap A_n \right | $$ The sets that are intersected are generated by SQL queries, and in an effort to keep things fast, I get a count of each query ahead of time, then take the set with the lowest count ($A_0$) and use those IDs as bounds on the rest of the big queries, so the intersection effectively becomes: $$ z = \left |\left ( A_0 \cap A_1 \right ) \cap \ldots \cap \left ( A_0 \cap A_n \right ) \right | $$ Even this strategy leaves me with some pretty big queries to run, since $\left | A_0 \right |$ can sometimes be large. My idea for dealing with that is taking a random sample of $A_0$ and intersecting it with the rest of the sets before extrapolating back to a proper estimate of $z$. **My question is: what is the best way to go about sampling and then extrapolating to get back to a value of $z$ that is, if not entirely accurate, has a predictable error range?** * * * Here's what I've tried so far (in pseudocode, sort of):               sample_threshold := 10000     factor := 1     if (len(A0) > sample_treshold) {         factor = sample_threshold / len(A0)     }          // Take a random sample of size 10000 from A0          // Intersect all the other sets with the A0 sample, then with each other     working_set := A0     for i, a := range A {         a = intersect(A0, a)         working_set = intersect(working_set, a)     }          z := len(working_set) * (1 / factor)      This code works, but seems to consistently overestimate `z`, with a lower sample size yielding a higher estimate. Additionally, I'm not sure how this would scale with more than two sets to intersect. I hope this question makes sense, let me know if I can clarify anything further. Also, if this question is off topic or belongs somewhere else, please let me know and I'm happy to move it. * * * Per Bill's comment, I ran some quick trials to show sample size vs. error. Each sample size bucket was run 20 times, and as you can see there's a pretty clear trend: ![Plot](http://i.stack.imgur.com/p5uc2.png)


title: Interpretation of interactions between categorical data

text: I have data that looks a bit like this                 A<-c(90,95,70,45,20,15)       B<-c(10,5,30,55,80,85)       X<-c("X1","X2","X1","X2","X1","X2")       Z<-c("Z1","Z1","Z2","Z2","Z3","Z3")      I've carried out a glm                mod1<-glm(cbind(A,B)~X*Z, binomial)     Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)   2.1972     0.3333   6.592 4.35e-11 ***     XX2           0.7472     0.5671   1.318 0.187659         ZZ2          -1.3499     0.3984  -3.388 0.000703 ***     ZZ3          -3.5835     0.4167  -8.600  < 2e-16 ***     XX2:ZZ2      -1.7952     0.6400  -2.805 0.005035 **      XX2:ZZ3      -1.0955     0.6801  -1.611 0.107230          The reference level is Z1 and X1 and I've interpreted the data in this way - As a main effect x is not significant but z is, with z2 and z3 significantly different form z1. We cannot tell the significance of the difference between z2 and z3 but we can work out the difference of the estimate... correct? However the x2:z2 interaction indicates x is having a different effect in z1 and z2? And the nonsignificant interaction of x2:z3 indicates that x is having the same effect in z1 and z3? Is this interpretation correct? Extra questions... How can I test the significance of the difference between z2 and z3, also the significance of the interaction between x2:z2 and x2:z3? Last one ... If the interaction is significant how do I interpret the significance of the main effect of say x?


title: How to calculate correlation between multiple?

text: I do have some 10-12 independent variable, I want to know how these 10-12 variables are related to a single variable(X). so I can decide that some variables are positively correlated to X and some are negatively correlated. I thought of doing correlation between X and one of the independent variables individually and decide. is there is a way to calculate multiple correlation.


title: Factor extracted more than expected

text: Can any one help me? I develop a scale with 31 items, and I expect to have five (5) factors in EFA but after conducting EFA, nine (9) factors were extracted. I want to know what cause this situation and how to solve the problem to have my 5 factors in my original model. Thank you.


title: Where does the offset go in Poisson/negative binomial regression?

text: (First of all, just to confirm, an offset variable functions basically the same way in Poisson and negative binomial regression, right?) Reading about the use of an offset variable, it seems to me that most sources recommend including that variable as an option in statistical packages (exp() in Stata or offset() in R). Is that functionally the same as converting your outcome variable to a proportion if you're modeling count data and there is a finite number that the count could have happened? My example is looking at employee dismissal, and I believe the offset here would simply be log(number of employees). And as an added question, I am having trouble conceptualizing what the difference is between these first two options (including exposure as an option in the software and converting the DV to a proportion) and including the exposure on the RHS as a control. Any help here would be appreciated.


title: Plotting a comet like animation for multiple variables

text: I am trying to code a visualization with 4 variables ( Carbon emission, Energy consumption, population and year) The data set i have collected so far looks like this With C1990 representing Carbon emission data 1990 and so on. ![enter image description here](http://i.stack.imgur.com/oMtR4.jpg) I have seen an online visualization whereby he uses some sort of animation which i suspect is comet3. His visualization looks like something below that i have recreated using paint. The animation will change in accordance to the year, the circle representing the countries and the size of the circle represent the carbon emission. Basically the circle will move and scale accordingly from year 1990 to 2010. I have been checking up on online resources but was unable to find something similar with a snipe code as sample. I am pretty new to matlab thus i am wondering if anyone here is able to get me off with a sample or guide me along ? ![enter image description here](http://i.stack.imgur.com/YCuxK.png)


title: How can I find the best linear combination of a set of matrices to approximate a target matrix?

text: I want to find $\theta$ such that $ \theta = argmin_{\theta} \left( \left|\left| Y - \sum_{i=1}^k \theta_i X_i \right|\right| \right) $ where $X_i$ and $Y$ are N x N matrices and $\theta$ is a weight vector that specifies how to linearly combine the $k$ $X$'s to approximate Y. This smells like a linear optimization problem though I'm unskilled in this kind of math and can't seem to formulate it as a linear program. I'd also be curious if anyone has any advice on how to learn this kind of math and/or problem formulation. Thanks!


title: Correlation with multiple x values for the same Y values

text: I had a quick question. I am trying to test the correlation between the % coral cover at each of my study sites, and the % of coral in sediment collected from those sites. The problem is that the coral cover per site is obtained independently of how the % coral in sediment is obtained. Specifically, I have a value for coral cover based on 3 independent samplings (linear transects) for each site. So I have the average coral cover per site. Then I collected samples of sand (3 at each site). For each of those samples, I get a % of coral in sediment value. So I wanted to see if per site, the amount of coral found in the sediment was related to the overall coral cover of the reef. Originally I had the % of sediment per sample as my X variables and the average coral cover per site as my Y variable. So ![enter image description here](http://i.stack.imgur.com/hIWyF.jpg) Notice though that I have the same Y value for each of the different x values per site. I wasn't sure if it was ok to do this. Someone told me that I needed to average the x values per site as well. When I did that, it leaves me with a sample size of only 6 (since I only have 6 reef sites to compare) and I am not able to identify any relationship. Does this make sense? What is the best way with this data to understand if there is any relationship. Thank you so much.


title: Geometric circle fitting with known radius

text: I have data points from a half circle and I already know the radius. I want to find the circle which best fits the points using a fixed radius. How can I do this? If I solve the problem using a typical circle fit algorithm the radius is too unstable due to "noise".


title: Converting SPSS syntax for regression to Stata commands

text: I am running a series of regressions to eliminate the effects of IQ, language ability and age from my variables (task performance). To do this, the variable is entered as the dependent variable (e.g. EFT in syntax below) in a regression with IQ, language ability and age as predictor variables, selecting for control group only (GROUP1 EQ 2 is the control group). This regression equation is then applied to the other groups (GROUP1 = 0 & 1) and the residuals are collected. The main problem I am having is converting the SPSS syntax below into Stata command (as I am relatively new to Stata). I need to use Stata to run the regression as my data are from related individuals and so I need to use the cluster command which isn't available in SPSS.               REGRESSION      /SELECT=GROUP1 EQ 2      /MISSING LISTWISE      /STATISTICS COEFF OUTS R ANOVA      /CRITERIA=PIN(.05) POUT(.10)     /NOORIGIN      /DEPENDENT EFT     /METHOD=ENTER AGE IQ LANG      /SAVE RESID (rEFT).


title: Logistic regression classifier with non-negative weights constraint

text: My feature data is defined in such a way that I believe all weights must be non-negative. I am looking for a reference discussing how to optimize the weights of a logistic regression classifier with the constraint that the weights must be non-negative, and perhaps also a constraint on the sign of the bias. Will replacing weights with squared weights and using the regular ML cost function with a local optimization scheme work?


title: When submitting analyses to FDA how does using R affect issues of software validation?

text: I have worked on clinical trials in FDA submissions for many years. I use SAS almost exclusively. Recently I discussed a consulting job I had to bootstrap a Deming regression. Bootstrapping is much more easily done in R I think and several experts have suggested ways that I can do this in R. Since this is going into a resubmission to the FDA for my client I am concerned as to whether software validation might become an issue. I hear that R is gaining acceptance with the FDA and experts tell me that traceability which is an attractive feature with SAS can be dealt with easily in R as well. My client did the original Deming regression with a package that I am not familiar with and validation was not an issue with the original submission. I want to help without adding any new headaches for them. Can the R experts here give me some guidance? If it will help I can mention some of the programming options that were suggested?


title: Representation of standard deviation in statistical range

text: I have computed the mean and standard deviation of a variable in a dataset, and I wish to represent these values in a report. Is the following the correct way to go about this? $1.23 \pm 0.52 \sigma$, where the mean is 1.23 and the standard deviation is 0.52. Alternatively, would a standard error be a more appropriate thing to report? The results are classification accuracies and the time taken to compute them.


title: Latest article or new development in cross validation?

text: Could anybody provide me a latest material related to Cross validation especially R package?


title: Error propagation SD vs SE

text: I have 3 to 5 measures of a trait per individual in two different conditions (A and B). I'm plotting the average for each individual in each condition and I use the standard error ( _i.e._ , $SD/\sqrt{N}$, with $N$ = number of measurements) as error bars. Now I want to plot the difference between the average measure per individual in condition A and condition B. I know I can determine the propagated error doing: $$SD=\sqrt{SD_A^2+SD_B^2}$$ but how can I propagate standard errors (since I'm dealing with averages of measurements) instead of standard deviations? Does this make sense at all?


title: Model selection with different outcome spaces

text: My question is about model selection/comparison in the case of discrete outcome spaces _and_ when the number of distinct outcomes depends on the dimension. (For the curious minded, this arises naturally, and is in fact the definition of dimension, in quantum mechanical models.) I'll give a simple analogy. Suppose I have a $d$-sided die I wish to know the bias of. But, moreover, I do not know $d$. If I've have seen $k$ distinct outcomes, then I know there is at least $k$ sides. But, perhaps there is more? My initial thought was to use model selection criteria, such as AIC. I could, of course, use AIC, but first I would have assume there is some maximum value of $k$. But this won't be practically relevant since the maximum likelihood term in the AIC will be the same regardless of how many additional sides I add to the die because the maximum likelihood estimator of the probability of those sides will be 0. Then, the model with the best AIC will be the one with the number of distinct observed sides: $k$. Seems overconfident to me. Moreover, I'd like a technique that doesn't require an assumption on the maximum number of possible outcomes. **My question, then** , is: are there model selection techniques other than AIC which can hedge between the absolute minimal description and ones that favor models which predict outcomes not yet observed?


title: Back-propagation in Neural Nets with >2 hidden layers

text: We have the following update formulas: Output layer (indexed by $k \in \{1, \dots, \text{number of classes} \}$). $o_k$ - value of output neuron $k$ $d_k$ - desired value of output neuron $k$ $\eta$ - learning rate $x_j$ - value of neuron on the last hidden layer $$\delta_k = o_k(1 - o_k)(d_k - o_k)$$ $$\omega_{jk} = \omega_{jk} + \Delta\omega_{jk}, \text{where } \Delta\omega_{jk} = \eta\delta_kx_j$$ Last hidden layer (HL) (indexed by $j \in \{1, \dots, \text{number of hidden neurons in the last HL} \}$): $$\delta_j = x_j(1 - x_j)\sum_k \omega_{jk}\delta_k$$ $$\omega_{ij} = \omega_{ij} + \Delta\omega_{ij}, \text{where } \Delta\omega_{ij} = \eta\delta_j x_i$$ **Does this update formulas hold generally for any number of hidden layers? I.e.:** $l$-th hidden layer ($l \in [0..L]$, where 0 - input layer, $L$ - output layer) (indexed by $p \in \{1, \dots, \text{number of hidden neurons in layer } l\},$ $j \in \{\text{ 1,..,number of hidden neurons in layer } l+1\}$ $$\delta_p = x_p(1 - x_p)\sum_j \omega_{pj}\delta_j$$ $$\omega_{pj} = \omega_{pj} + \Delta\omega_{pj}, \text{where } \Delta\omega_{pj} = \eta\delta_j x_p$$ If it is true, it's possible to recursively compute all $\delta$s for all neurons (except the input ones) starting with output and going deeper to the hidden nodes and then simply modify the each weight based on found $\delta$ and $x$. Correct?


title: Reordering a dendrogram in R

text: I have a dendrogram called `x` and I want to reorder it. First, I make sure that `x` is a dendrogram, and it is:               > class(x)      [1] "dendrogram"      Then I check the order of things in my dendrogram, and it's exactly as it is drawn, i.e. everything makes sense:               > order.dendrogram(x)      [1] 6 7 8 5 9 10 11 1 3 2 4 13 16 12 14 15 17      Now I decide on a new order:               > y <- c(6,7,8,5,9,1,3,2,4,13,16,17,15,14,12,10,11)      I now try to reorder my dendrogram and store the reordered one as `rx`:               > rx <- reorder(x, y, agglo.FUN=mean)      I check the order of `rx`, and it is reordered from `x`, but it is not reordered in quite the way I want (almost, but not quite):               > order.dendrogram(rx)      [1] 6 7 8 9 5 4 2 1 3 16 13 17 15 14 12 10 11      I played around a bit with the `agglo.FUN=mean` argument, but never could get the dendrogram to be ordered just like I wanted. The ordering I'm asking for is completely compatible with the dendrogram - I'm just rotating some branches. Does anyone know what I'm doing wrong? Thanks. Eric


title: Make R report error on using non-existent column name in a data frame

text: In our usage of R for a non-trivial data analysis and estimation project, we've been repeatedly burnt by how tolerant R is toward misspelled or missing columns in a data frame. Typical example is calculating the weighted mean of a variable MYVAR in a data frame using another variable WEIGHT for weights:               m <- weighted.mean(tbl$MYVAR, w = tbl$WEIGHT, na.rm = TRUE)      Suppose I make a typo in the WEIGHT name in the operation above. What will happen in that R will expand my misspelled column into NULL and will use it for performing the weighted mean resulting in a non-weighted one. Therefore, the question: is there any way to make R treat attempts to "read" a non-existent variables in a data frame as an error?


title: Estimating the probability of a person getting a question right

text: We need to estimate the probability of a person getting a question right for a given content area given his history of getting questions in the same content area right in the past. We would also presumably have records on how others have done on this question and in this content area. Is there a good way or ways of doing this? (This is sort of an education theory question, but I couldn't find a better place to post this question.)


title: Why is pruning not needed for random forest trees?

text: Breiman says that the trees are grown with out pruning. Why? I mean to say that there must be a solid reason why the trees in random forest are not pruned. On the other hand it is considered very important to prune a single decision tree to avoid over fitting. Is there some literature available to read for this reason? Of course the trees may not be correlated but there still will be a possibility of over fitting.


title: Comparing  overdispersion distributions

text: I am dealing with handling overdispersed count data(Poisson distributions fails to fit). I need to compare three different mixtures (Gamma, Log-Normal and Inverse Gaussian) of Poisson rate parameter distributions using a simulation study. I would like to know the available methods(apart from Likelihood ration and Pearson's Chi Square test) to compare three completely different distributions to fit overdispersed count data. Any reference also highly appreciated!! Thank you for your kindness.


title: Chi-squared and repeated measures

text: I would like to determine whether there is a deviation in the male:female ratio of a certain strain of laboratory mouse from the normal situation (with male:female = 1:1). I have several mouse pairs, and several litters for each pair. So, for instance:               Pair 1        Litter 1  5M:3F        Litter 2  3M:4F        Litter 3  7M:4F        ...     Pair 2        Litter 1  4M:5F        Litter 2  2M:5F        Litter 3  6M:2F     ...      What would be the best way to analyse this situation?   I thought of using chi-squared, but I am not sure about how to handle the fact that I have repeated measures for each pair.


title: Variance of two weighted random variables

text: Let: Standard deviation of random variable $A =\sigma_{1}=5$ Standard deviation of random variable $B=\sigma_{2}=4$ Then the variance of A+B is: $Var(w_{1}A+w_{2}B)= w_{1}^{2}\sigma_{1}^{2}+w_{2}^{2}\sigma_{2}^{2} +2w_{1}w_{2}p_{1,2}\sigma_{1}\sigma_{2}$ Where: $p_{1,2}$ is the correlation between the two random variables. $w_{1}$ is the weight of random variable A $w_{2}$ is the weight of random variable B $w_{1}+w_{2}=1$ The figure below plots the variance of A and B as the weight of A changes from 0 to 1, for the correlations -1 (yellow),0 (blue) and 1 (red). ![alt text](http://i.stack.imgur.com/ATJSF.png) How did the formula result in a straight line (red) when the correlation is 1? As far as I can tell, when $p_{1,2}=1$, the formula simplifies to: $Var(w_{1}A+w_{2}B)= w_{1}^{2}\sigma_{1}^{2}+w_{2}^{2}\sigma_{2}^{2} +2w_{1}w_{2}\sigma_{1}\sigma_{2}$ How can I express that in the form of $y=mx+c$? Thank you.


title: Summary and frequency tables in JMP

text: 1. Does anyone know how to create a nice looking summary table in JMP for several variables? I use Tables>Summary and the table I get has one row that includes for example mean(var1) mean(var2) mean(var3) median(var1) median(var2) median(var3) and so on. I would like each variable to have row of its own, with columns representing the statistics.   2. Does anyone know how to create a frequency table when creating a histogram in JMP? I get a histogram by using Analyze>Distributions but I don't see any option of getting the frequency table. Thank you, Naama Wolf


title: Expectation of product of Gaussian random variables

text: Say we have two Gaussian random vectors $p(x_1) = N(0,\Sigma_1), p(x_2) = N(0,\Sigma_2)$, is there a well known result for the expectation of their product $E[x_1x_2^T]$ without assuming independence?


title: Confidence interval about a multivariate normal probability

text: From a given sample $x_i \underset{\mathrm{iid}}{\sim} {\cal N}(\mu, \sigma^2)$ it is possible to get a confidence interval about the probability $\Pr(x_i \geq a)$ for any number $a$, by "inverting" some tolerance intervals. If ${\boldsymbol x}_i \underset{\mathrm{iid}}{\sim}{\cal N}_p({\boldsymbol \mu}, \Sigma)$ is a sample from a $p$-variate normal distribution, how to get a confidence interval about the probability $\Pr(x_{i1} \geq a_1, \ldots, x_{ip} \geq a_p)$? It is easy to get a credible interval in the Bayesian framework, but what about frequentist methods? Is the bootstrap asymptotically valid? Does there exist another known approach? **EDIT 15/09/2012:** The Delta method should be possible. If someone is comfortable with the calculations I would be glad if he posts the solution here ;)


title: How to get unbiased estimate for error variance in ordinary linear model

text: I would like to estimate the variance of the error term in the ordinary linear regression model. The obvious estimate is the sample variance of the residuals, however this estimate consistently underestimates the real value. A simple simulation in R confirms that:                   x = rep(0:1, c(10, 10))         a = replicate(1000, var(residuals(lm(rnorm(20, mean = x, sd = 1) ~ x)))); mean(a)         > 0.9418632      This is of course makes sense, since when optimizing a linear model we are actually optimizing this error and the real generating model would always give somewhat bigger errors. However, is there a good way to get an unbiased estimate for this error? One possibility is probably using some kind of cross- validation scheme, but is there any way to calculate this analytically?